# C-3PR
# ========
# ### Evaluating Scientific **C**laims: **P**ost **P**ublication **P**eer **R**eview, **R**eplication and **R**eproduction
#
# An R function library created in the contex of analysing and visualising data generated by several Open Science intiatives.
#
# *C-3PR Created by [Fred Hasselman](http://fredhasselman.com) on behalf of the ManyLabs1, ManyLabs2, RPP and Curate Science projects*
#
# Some functions included are based on code by others (adjusted in varying degrees of mutilation)

# SOURCE FROM GITHUB -------------------------------------------------
#
# Use this code (devtools package) to source it directly from GitHub:
# require(devtools)
# source_url("https://raw.githubusercontent.com/FredHasselman/toolboxR/master/C-3PR.R")

# Package install / load / unload -----------------------------------------

# Packages in the list argument need will be installed if necessary and loaded
in.IT <- function(need=NULL,inT=TRUE){
    ip <- .packages(all.available=TRUE)
    if(any((need %in% ip)==FALSE)){
        if(inT==TRUE){
            install.packages(need[!(need %in% ip)])
        } else {
            cat('Package(s):\n',paste(need[(need %in% ip)==FALSE],sep='\n'),'\nnot installed.\nUse in.IT(c("packagename1","packagename2",...),inT=TRUE)')
            need <- need[(need %in% ip)==TRUE]
        }
    }
    ok <- sapply(1:length(need),function(p) require(need[[p]],character.only=TRUE))
}

# Packages in the list argument loose will be unloaded, if necessary (unT=T) uninstalled

un.IT <- function(loose,unT=FALSE){
    dp <- .packages()
    if(any(loose %in% dp)){
        for(looseLib in loose[(loose %in% dp)]){detach(paste0("package:",looseLib), unload=TRUE,character.only=TRUE)}
    }
    rm(dp)
    if(unT==TRUE){
        dp <- .packages(all.available=TRUE)
        if(any(loose %in% dp)){remove.packages(loose[(loose %in% dp)])}
    }
}

in.IO <- function(){
    # I/O and data handling tools
    in.IT(c("xlsx","plyr","doBy","reshape2","RCurl","XML","httr","dplyr"))
}

in.PAR <- function(){
    # Parallel computing tools
    in.IT(c("parallel","doParallel","foreach"))
}

in.PLOT <- function(useArial = F,afmPATH="~/Dropbox/Public"){
    # Load packages for plotting with default option to setup Arial as the pdf font for use in figures.
    if(useArial==T){
        # Set up PDF device on MAC OSX to use Arial as a font in Graphs
        set.Arial(afmPATH)
    }
    in.IT(c("lattice","latticeExtra","gplots","ggplot2","grid","gridExtra","scales","beanplot","effects","RColorBrewer"))
}

# DATA loading and cleaning -----------------------------------------------
df.Clean <- function(df,Sep="."){
    in.IT('dplyr')
    nms   <- colnames(df)
    rws   <- rownames(df)

    # Change punctuation and blankss in variable names to points
    nmsP  <- gsub("([[:punct:]]|[[:blank:]])+","+",nms)
    nmsPP <- gsub("(^[+]|[+]$)+","",nmsP)
    nmsPP <- gsub("[+]",Sep,nmsPP)
    # Check for double names
    ifelse(length(unique(nmsPP))==length(nmsPP),{nms <- nmsPP},{
        id2 <- which(plyr::laply(nmsPP,function(n) sum(nmsPP%in%n))>1)
        nms <- nmsPP
        nms[id2] <- paste(nmsPP[id2],id2,sep=".")})

    colnames(df) <- nms
    df      <- dplyr::select(df,which(nms%in%nms[nms!=""]))
    df[ ,1] <- paste0("Row.",seq(1,nrow(df)))
    colnames(df)[1] <- paste("Local","ID",sep=Sep)
    return(list(df=df,
                nms=nms,
                rws=rws))
}

get.GoogleSheet <- function(url=NULL,data=c('ML1data','ML2masteRkey','RPPdata')[3],dfCln=FALSE,Sep = "."){
    in.IT(c('dplyr','httr'))
    if(is.null(url)){
        switch(data,
               ML1data        = url <- 'https://docs.google.com/spreadsheets/d/19ay71M8jiqIZhSj3HR0vqaJUwAae1QHzBybjBu5yBg8/export?format=csv',
               ML2masteRkey   = url <- 'https://docs.google.com/spreadsheets/d/1fqK3WHwFPMIjNVVvmxpMEjzUETftq_DmP5LzEhXxUHA/export?format=csv',
               RPPdata        = url <- ''
        )}
    # GET(url) will only get 100 rows, thanks to Sacha Epskamp for this "complete scrape" code.
    tmp  <- tempfile()
    info <- GET(url, write_disk(tmp, overwrite = TRUE))
    df   <- tbl_df(read.csv(tmp,stringsAsFactors=F, fill = T,header=T))
    if(dfCln==TRUE){df   <- df.Clean(df)} else {df$df <- df}

    return(list(df = df$df,
                info = list(Info=info,
                            GoogleSheet.colnames=tbl_df(data.frame(ori.colnames=df$nms)),
                            GoogleSheet.rownames=tbl_df(data.frame(ori.rownames=df$rws)))))
}


get.OSFfile <- function(# Function to download OSF file modified from code by Sacha Epskamp
    code,  #Either "https://osf.io/XXXXX/" or just the code
    dir = tempdir(), # Output location
    scanMethod, #  "readLines" or "RCurl". Leave missing to automatically chose
    downloadMethod = c("httr","downloader","curl"), # First one is chosen
    dataFrame = TRUE,
    sep = ',',
    dfCln = FALSE
){
    # Check if input is code:
    if (!grepl("osf\\.io",code)){
        URL <- sprintf("https://osf.io/%s/",code)
    } else URL <- code

    # Scan page:
    if (grepl("Windows",Sys.info()[['sysname']],ignore.case=TRUE)){
        try(setInternet2(TRUE))
    }

    if (missing(scanMethod)){
        scanMethod <- ifelse(grepl("Windows",Sys.info()[['sysname']],ignore.case=TRUE), "readLines", "RCurl")
    }
    if (scanMethod == "readLines"){
        Page <- paste(readLines(URL),collapse="\n")
    } else if (scanMethod == "RCurl"){
        library("RCurl")
        Page <- RCurl::getURL(URL)
    } else if (scanMethod == "httr"){
        Page <- httr::GET(URL)
        Page <- paste(Page,collapse="\n")
    } else stop("Invalid scanMethod")

    # Extract download link(s):
    #   Link <- regmatches(Page, regexpr("(?<=download: \\').*?\\?action=download(?=\\')", Page, perl = TRUE))
    #
    #   # Stop if no download link:
    #   if (length(Link)==0){
    #     stop("No download link found")
    #   }
    #   # (just in case) if more than one, warning:
    #   if (length(Link)>1){
    #     warning("Multiple download links found, only first is used")
    #     Link <- Link[1]
    #   }
    # Create download link:
    URL <- gsub("/$","",URL)
    #   Link <- paste0(URL,"/?action=download&version=1")
    Link <- paste0(URL,"/?action=download")

    # Extract file name:
    FileName <- regmatches(Page,gregexpr("(?<=\\<title\\>OSF \\| ).*?(?=\\</title\\>)", Page, perl=TRUE))[[1]]
    FullPath <- paste0(dir,"/",FileName)

    info <- NULL
    # Download file:
    if (downloadMethod[[1]]=="httr"){
        library("httr")
        info <- httr::GET(Link, httr::write_disk(FullPath, overwrite = TRUE))
    } else if (downloadMethod[[1]]=="downloader"){
        library("downloader")
        downloader:::download(Link, destfile = FullPath, quiet=TRUE)
    } else if (downloadMethod[[1]]=="curl"){
        system(sprintf("curl -J -L %s > %s", Link, FullPath), ignore.stderr = TRUE)
    }  else stop("invalid downloadMethod")

    df <- NULL
    if(dataFrame==TRUE){
        if(grepl('xls',FileName)){
            df <- tbl_df(read.xlsx2(file=FullPath,sheetIndex=1))
        } else {
            df <- tbl_df(read.table(FullPath,stringsAsFactors=F,fill = T,header=T,sep=sep, comment.char = "",quote = "\""))
        }
        if(dfCln==TRUE){df <- df.Clean(df)} else {df$df <- df}

        return(list(df   = df$df,
                    info = list(FilePath=FullPath,
                                Info=info,
                                ori.Colnames=tbl_df(data.frame(ori.colnames=df$nms)),
                                ori.Rownames=tbl_df(data.frame(ori.rownames=df$rws))
                    )))
    } else {
        # Return location of file:
        return(FilePath=FullPath)
    }
}


get.Order <- function(df, S1=TRUE){

    ifelse(S1,{
        url <-  "https://docs.google.com/spreadsheets/d/1al8b5nv9AoNdOOlI-RfXf5bCTJY2ophdUIiYhCFGegI/pub?gid=269287357&single=true&output=csv"
        cols <- 3:15
    },{
        url <- "https://docs.google.com/spreadsheets/d/1AvJguRhN8i7MsbcjnGtEnNZ9b-LNTNKINSoHm8Z2f28/pub?gid=1370595041&single=true&output=csv"
        cols <- 3:17
        }
    )

    df.Order <- get.GoogleSheet(url = url)$df

    ProblemID <- list()
    cnt = 0

    for(i in 1:nrow(df)){

        if((nchar(df$StudyOrder[i])==0)|(df$IDiffOrder[i]=="")){
            cnt = cnt + 1
            if((nchar(df$StudyOrder[i])==0)&(df$IDiffOrder[i]=="")){
                df$StudyOrderN[i] <- NA
                df$IDiffOrderN[i] <- NA
                Problem = paste("No study order strings in",df$StudyOrder[i],"and",df$IDiffOrder[i])

            } else {

                if(nchar(df$StudyOrder[i])==0){ df$StudyOrderN[i] <- NA}
                if(df$IDiffOrder[i]==""){
                    df$IDiffOrderN[i] <- NA
                    Problem = paste("Wrong study order var? Ind.Diff var",df$IDiffOrder[i],"is empty")
                }
            }

            ProblemID[[cnt]] <- cbind(rowNum     = i,
                                      fileName   = df$.id[i],
                                      ResponseID = df$ResponseID[i],
                                      Problem    = Problem,
                                      StudyOrder = df$StudyOrder[i],
                                      IDiffOrder = df$IDiffOrder[i])
        } else {

            OrderString <- unlist(strsplit(x = df$StudyOrder[i], split = "[|]"))
            ls   <- list()
            Problem.s <- list()
            cnt2 <- 0

            for(s in 1:length(OrderString)){
                ls[[s]] <- colnames(df.Order)[cols][df.Order[df.Order$Filename%in%df$.id[i], cols] %in% OrderString[[s]]]
                if(length(ls[[s]])==0){
                    cnt2 = cnt2 + 1
                    if(OrderString[[s]]%in%df.Order[df.Order$Filename%in%df$.id[i], cols]){
                        Problem.s[[cnt2]] <- paste(OrderString[[s]],"(", colnames(df.Order)[cols][df.Order[df.Order$Filename%in%df$.id[i], cols] == OrderString[[s]]], ")")
                        ls[[s]] <- NA
                    } else {
                        Problem.s[[cnt2]] <- paste(OrderString[[s]], "( mismatch )")
                    }
                }
            }

            if(cnt2!=0){
                cnt = cnt + 1
                ProblemID[[cnt]] <- cbind(rowNum     = i, fileName = df$.id[i],
                                          ResponseID = df$ResponseID[i],
                                          Problem    = paste(unlist(Problem.s), collapse="|"),
                                          StudyOrder = paste(unlist(ls) ,collapse="|"),
                                          IDiffOrder = df$IDiffOrder[i])
            }

            df$StudyOrderN[i] <- paste(unlist(ls) ,collapse="|")

        }

        df$IDiffOrderN[i] <- df$IDiffOrder[i]

    }

    return(list(df = df,
                Problems = ldply(ProblemID))
           )

}


# PLOTS -------------------------------------------------------------------
disp <- function(message='Hello world!', header = TRUE, footer = TRUE){

    mWidth <- max(laply(message,nchar))

    if(is.character(header)){
        hWidth <- max(laply(header,nchar))
        mWidth <- max(hWidth,mWidth)
    }

    dmessage <- list()
    for(m in 1:length(message)){
       # b <- floor((mWidth-nchar(message[m]))/2)
        e <- mWidth-nchar(message[m])
        dmessage[[m]] <- paste0('§ ',message[m]) #,paste0(rep(' ',e),collapse=""),'\n\t')
                                #paste0('§ ',paste0(rep(" ",mWidth),collapse=""),' §'))
    }
    # if(m > 1){dmessage[[m]] <- paste0(dmessage[[m]],}

   # mWidth <- max(laply(dmessage, nchar))
    banner <- paste0(rep('~', mWidth), collapse = "")
    if(is.character(header)){
        b <- floor((nchar(banner)-nchar(header))/2)
        e <- ceiling((nchar(banner)-nchar(header))/2)
            leader <- paste0('\n\t',paste0(rep('~',b),collapse=""),header,paste0(rep('~',e),collapse=""))
        }
    if(header == TRUE){
            leader <- banner
        }
    if(header == FALSE){
            leader <- paste0('§') #,paste0(rep(" ",nchar(banner)-2),collapse="")) #,'§')
        }

    if(footer){
            cat(paste0('\n\t',leader,'\n\t',dmessage,'\n\t',banner,'\n'))
        } else {
            cat(paste0('\n\t',leader,'\n\t',dmessage))
        }
}


gg.theme <- function(type=c("clean","noax")[1],useArial = F, afmPATH="~/Dropbox"){
    require(ggplot2)
    if(useArial){
        set.Arial(afmPATH)
        bf_font="Arial"
    } else {bf_font="Helvetica"}

    switch(type,
           clean = theme_bw(base_size = 16, base_family=bf_font) +
               theme(axis.text.x     = element_text(size = 14),
                     axis.title.y    = element_text(vjust = +1.5),
                     panel.grid.major  = element_blank(),
                     panel.grid.minor  = element_blank(),
                     legend.background = element_blank(),
                     legend.key = element_blank(),
                     panel.border = element_blank(),
                     panel.background = element_blank(),
                     axis.line  = element_line(colour = "black")),

           noax = theme(line = element_blank(),
                        text  = element_blank(),
                        title = element_blank(),
                        plot.background = element_blank(),
                        panel.border = element_blank(),
                        panel.background = element_blank())
    )
}

plotHolder <- function(useArial = F,afmPATH="~/Dropbox"){
    require(ggplot2)
    ggplot() +
        geom_blank(aes(1,1)) +
        theme(line = element_blank(),
              text  = element_blank(),
              title = element_blank(),
              plot.background = element_blank(),
              #           panel.grid.major = element_blank(),
              #           panel.grid.minor = element_blank(),
              panel.border = element_blank(),
              panel.background = element_blank()
              #           axis.title.x = element_blank(),
              #           axis.title.y = element_blank(),
              #           axis.text.x = element_blank(),
              #           axis.text.y = element_blank(),
              #           axis.ticks = element_blank()
        )
}

set.Arial <- function(afmPATH="~/Dropbox"){
    # Set up PDF device on MAC OSX to use Arial as a font in Graphs
    if(nchar(afmPATH>0)){
        if(file.exists(paste0(afmPATH,"/Arial.afm"))){
            Arial <- Type1Font("Arial",
                               c(paste(afmPATH,"/Arial.afm",sep=""),
                                 paste(afmPATH,"/Arial Bold.afm",sep=""),
                                 paste(afmPATH,"/Arial Italic.afm",sep=""),
                                 paste(afmPATH,"/Arial Bold Italic.afm",sep="")))
            if(!"Arial" %in% names(pdfFonts())){pdfFonts(Arial=Arial)}
            if(!"Arial" %in% names(postscriptFonts())){postscriptFonts(Arial=Arial)}
            return()
        } else {disp(header='useArial=TRUE',message='The directory did not contain the *.afm version of the Arial font family')}
    } else {disp(header='useArial=TRUE',message='Please provide the path to the *.afm version of the Arial font family')}
}

#' Graphical display of a textual table.
#'
#' @param d data.frame or matrix
#' @param widths optional vector to specify column widths
#' @param heights optional vector to specify row heights
#' @param fg.par control parameters for text grobs
#' @param bg.par control parameters for rect grobs
#' @param padding unit of length 2
#' @export
#' @examples
#' \donttest{
#' d <- head(iris, 3)
#' core <- gtable_table(d,
#'                      fg.par = list(col=1:5, fontsize=c(10,12,15)),
#'                      bg.par = list(fill=1:2, alpha=0.5))
#' colhead <- gtable_table(t(colnames(d)))
#' rowhead <- gtable_table(c("", rownames(d)))
#' g <- rbind(colhead, core)
#' g <- cbind(rowhead, g)
#' grid.newpage()
#' grid.draw(g)
#' g2 <- gtable_table(t(colnames(d)), widths = unit(rep(1, ncol(d)), "null"))
#' grid.newpage()
#' grid.draw(g)
#' }
gtable_table <- function(d, widths, heights,
                         fg.par = list(col = "black"),
                         bg.par = list(fill = NA),
                         padding = unit(c(4, 4), "mm")){

    label_matrix <- as.matrix(d)

    nc <- ncol(label_matrix)
    nr <- nrow(label_matrix)
    n <- nc*nr

    fg.par <- lapply(fg.par, rep, length.out = n)
    bg.par <- lapply(bg.par, rep, length.out = n)

    fg.param <- data.frame(fg.par, label = as.vector(label_matrix),
                           stringsAsFactors=FALSE)
    bg.param <- data.frame(bg.par, id = seq_len(n),
                           stringsAsFactors=FALSE)

    labels <- plyr::mlply(fg.param, cell_content)
    backgrounds <- plyr::mlply(bg.param, cell_background)

    label_grobs <- matrix(labels, ncol = nc)

    ## some calculations of cell sizes
    row_heights <- function(m){
        do.call(unit.c, apply(m, 1, function(l)
            max(do.call(unit.c, lapply(l, grobHeight)))))
    }

    col_widths <- function(m){
        do.call(unit.c, apply(m, 2, function(l)
            max(do.call(unit.c, lapply(l, grobWidth)))))
    }

    if(missing(widths))
        widths <- col_widths(label_grobs) + padding[1]
    if(missing(heights))
        heights <- row_heights(label_grobs) + padding[2]

    ## place labels in a gtable
    g <- gtable_matrix("table", grobs=label_grobs,
                       widths = widths,
                       heights = heights)

    ## add the background
    g <- gtable_add_grob(g, backgrounds, t=rep(seq_len(nr), each=nc),
                         l=rep(seq_len(nc), nr), z=0, name="fill")

    g
}


cell_content <- function(...){
    dots <- list(...)
    gpar.names <- c("col", "cex", "fontsize", "lineheight",
                    "font", "fontfamily", "fontface", "alpha")
    other.names <- c("label", "hjust", "vjust", "rot", "x", "y")
    gpar.args <- dots[intersect(names(dots), gpar.names)]
    gp <- do.call(gpar, gpar.args)
    other.args <- dots[intersect(names(dots), other.names)]
    do.call(textGrob, c(other.args, list(gp = gp)))

}

cell_background <- function(...){

    dots <- list(...)
    gpar.names <- c("fill", "col", "lty", "lwd", "cex", "alpha",
                    "lineend", "linejoin", "linemitre",
                    "lex")
    gpar.args <- dots[intersect(names(dots), gpar.names)]
    gp <- do.call(gpar, gpar.args)
    do.call(rectGrob, list(gp = gp))

}


gtable_arrange <- function(..., grobs=list(), as.table=TRUE,
                           top = NULL, bottom = NULL,
                           left = NULL, right = NULL, draw=TRUE){
    require(gtable)
    # alias
    gtable_add_grobs <- gtable_add_grob

    dots <- list(...)
    params <- c("nrow", "ncol", "widths", "heights",
                "respect", "just", "z") # TODO currently ignored

    layout.call <- intersect(names(dots), params)
    params.layout <- dots[layout.call]

    if(is.null(names(dots)))
        not.grobnames <- FALSE else
            not.grobnames <- names(dots) %in% layout.call

    if(!length(grobs))
        grobs <- dots[! not.grobnames ]

    ## figure out the layout
    n <- length(grobs)
    nm <- n2mfrow(n)

    if(is.null(params.layout$nrow) & is.null(params.layout$ncol))
    {
        params.layout$nrow = nm[1]
        params.layout$ncol = nm[2]
    }
    if(is.null(params.layout$nrow))
        params.layout$nrow = ceiling(n/params.layout$ncol)
    if(is.null(params.layout$ncol))
        params.layout$ncol = ceiling(n/params.layout$nrow)

    if(is.null(params.layout$widths))
        params.layout$widths <- unit(rep(1, params.layout$ncol), "null")
    if(is.null(params.layout$heights))
        params.layout$heights <- unit(rep(1,params.layout$nrow), "null")

    positions <- expand.grid(row = seq_len(params.layout$nrow),
                             col = seq_len(params.layout$ncol))
    if(as.table) # fill table by rows
        positions <- positions[order(positions$row),]

    positions <- positions[seq_along(grobs), ] # n might be < ncol*nrow

    ## build the gtable, similar steps to gtable_matrix

    gt <- gtable(name="table")
    gt <- gtable_add_cols(gt, params.layout$widths)
    gt <- gtable_add_rows(gt, params.layout$heights)
    gt <- gtable_add_grobs(gt, grobs, t = positions$row,
                           l = positions$col)

    ## titles given as strings are converted to text grobs
    if (is.character(top))
        top <- textGrob(top)
    if (is.character(bottom))
        bottom <- textGrob(bottom)
    if (is.character(right))
        right <- textGrob(right, rot = -90)
    if (is.character(left))
        left <- textGrob(left, rot = 90)

    if(!is.null(top)){
        gt <- gtable_add_rows(gt, heights=grobHeight(top), 0)
        gt <- gtable_add_grobs(gt, top, t=1, l=1, r=ncol(gt))
    }
    if(!is.null(bottom)){
        gt <- gtable_add_rows(gt, heights=grobHeight(bottom), -1)
        gt <- gtable_add_grobs(gt, bottom, t=nrow(gt), l=1, r=ncol(gt))
    }
    if(!is.null(left)){
        gt <- gtable_add_cols(gt, widths=grobWidth(left), 0)
        gt <- gtable_add_grobs(gt, left, t=1, b=nrow(gt), l=1, r=1)
    }
    if(!is.null(right)){
        gt <- gtable_add_cols(gt, widths=grobWidth(right), -1)
        gt <- gtable_add_grobs(gt, right, t=1, b=nrow(gt), l=ncol(gt), r=ncol(gt))
    }

    if(draw){
        grid.newpage()
        grid.draw(gt)
    }
    invisible(gt)
}


#function to create geom_ploygon calls
fill_viol<-function(gr.df,gr,qtile,probs){
    # SETUP VIOLIN QUANTILE PLOTS -----------------------------------
    # This is adapted from: http://stackoverflow.com/questions/22278951/combining-violin-plot-with-box-plot
    #   g.df <- rbind(g.df[1,],g.df,g.df[nrow(g.df),])
    #   g.df$y[1]         <-min(qtile)
    #   g.df$y[nrow(g.df)]<-max(qtile)


    #fill_viol<-function(v,gr){
    #   quants<-mutate(v,
    #                  x.l=x-violinwidth/2,
    #                  x.r=x+violinwidth/2,
    #                  cuts=cut(y,quantile(df[df$grp==gr,"val"])))
    #   # add 1/2 width each way to each x value
    #   plotquants<-data.frame(x=c(quants$x.l,rev(quants$x.r)),   # left x bottom to top, then right x top to bottom
    #                          y=c(quants$y,rev(quants$y)),       # double up the y values to match
    #                          id=c(quants$cuts,rev(quants$cuts)))# cut by quantile to create polygon id
    #   geom_polygon(aes(x,y,fill=as.factor(id)),data=plotquants) # return the geom_ploygon object
    # }

    #cuts <- cut(g.df$y, breaks = breaks, include.lowest=T, right=T)
    #levels(cuts)
    ifelse(is.null(qtile),{
        cuts <- cut(gr.df$y, breaks = quantile(gr.df$y, probs, na.rm=T, type=3, include.lowest = T, right = T), na.rm=T)},{
            cuts <- cut(gr.df$y, breaks = qtile, na.rm=T)
        }
    )
    #qtile<-qtiles[gr,]

    quants <- mutate(gr.df,
                     x.l=x-violinwidth/2,
                     x.r=x+violinwidth/2,
                     cuts=cuts)

    plotquants <- data.frame(x=c(quants$x.l,rev(quants$x.r)),
                             y=c(quants$y,rev(quants$y)),
                             id=c(quants$cuts,rev(quants$cuts)))

    #cut by quantile to create polygon id
    geom <- geom_polygon(aes(x=x,y=y,fill=factor(id)),data=plotquants,alpha=1)

    return(list(quants=quants,plotquants=plotquants,geom=geom))
}

vioQtile <- function(gg=NULL,qtiles=NULL,probs=seq(0,1,.25),labels=paste(probs[-1]*100),withData=FALSE){
    require(ggplot2)
    # SETUP VIOLIN QUANTILE PLOTS -----------------------------------
    # This is adapted from: http://stackoverflow.com/questions/22278951/combining-violin-plot-with-box-plot
    #
    # Changed:
    # - Deal with 'empty' quantile groups
    # - Deal with original data
    # - More input, more output
    g.df <- ggplot_build(gg)$data[[1]]    # use ggbuild to get the outline co-ords

    #   gg <- gg +
    #     geom_dotplot(aes(y=EffectSize),fill="grey80", color="grey90", stackdir='center', binaxis='y', dotsize=.4, alpha=.9, stackratio=1.4, binpositions="all", method="histodot", na.rm=T, binwidth=.05)
    #

    ifelse(is.null(qtiles),{
        gg <- gg + lapply(unique(g.df$group), function(x) fill_viol(g.df[g.df$group==x, ],x,NULL,probs)$geom)},{
            gg <- gg + lapply(unique(g.df$group), function(x) fill_viol(g.df[g.df$group==x, ],x,qtiles[x, ],probs)$geom)}
    )

    #lapply(unique(g.df$group), function(x) fill_viol(g.df[g.df$group==x,],qtiles[x,],probs)$geom) +

    gg <- gg + geom_hline(aes(yintercept=0)) +
        scale_fill_grey(name="Quantile\n",labels=labels,guide=guide_legend(reverse=T,label.position="right")) +
        stat_summary(fun.y=median, geom="point", size=8, color="grey80",shape=21,fill="white")

    #     scale_fill_continuous(name="Quantile\n",labels=paste(probs[-1]*100),guide=guide_legend(reverse=T),low="grey20",high="grey80")
    #     scale_fill_gradient2(low = muted("sienna4",l=60),
    #                          mid = muted("thistle4",l=60),
    #                          high = muted("lightsteelblue4",l=60),
    #                          midpoint =0)

    if(withData){
        ifelse(is.null(qtiles),{
            ggData <- lapply(unique(g.df$group), function(x) fill_viol(g.df[g.df$group==x,],x,NULL,probs))},{
                ggData <- lapply(unique(g.df$group), function(x) fill_viol(g.df[g.df$group==x,],x,qtiles[x,],probs))
            }
        )
        return(list(ggGraph=gg,ggData=ggData))
    } else {
        return(gg)
    }
}



# MANYLABS 2 --------------------------------------------------------------
get.results <- function(key,dfr,s,group=NULL){

    dfr <- tbl_df(dfr)
    disp('Get info to create a dataset for the current study',header='get.results')
    inf <- get.info(key[s,],colnames(dfr))

    disp('# Generate chain to select variables for the data frame and create a filter chain for the variables to use for analysis\n')
    # Info based on KeyTable information in study.vars, cases.include, site.include, params.NA
    fltr <- get.chain(inf)

    cat('# Apply the df chain to select relevant subset of variables\n')
    dfr <- eval(parse(text=paste("dfr",fltr$df)))

    cat('# Get a list containing the data frames to be used in the analysis\n')
    ifelse(is.null(group),{
        sr <- get.sourceData(fltr,dfr,inf)
        sr <- list(sr)},{
            grp <- unique(eval(parse(text=paste0("dfr$",group))))
            sr <- llply(grp, function(sID) get.sourceData(fltr,filter(dfr,.id==sID),inf))
        })

    cat("# Organize and Calculate variables for the analysis using function according to ML2.info[s,'stat.vars']\n")
    var <- llply(sr, function(sID){
        test <- try.CATCH(eval(parse(text=paste0(key[s,'stat.vars'],'(sID)',collapse=""))))
        if(is.null(test$warning)){
            return(test$value)
        } else {
            cat(s,key[s,'study.name'],' stat.test failed\n')
        }
    })

    cat("# Run the analysis in 'stat.test': ",key[[s,'stat.test']],"\n")
    stat.params <- inf$stat.params
    stat.test   <- llply(var, function(vID){
        test <- try.CATCH(with(vID,eval(parse(text=key[s,'stat.test']))))
        if(is.null(test$warning)){
            return(test$value)
        } else {
            cat(s,key[[s,'study.name']],' stat.test failed\n')
        }
    })

    return(list(info=inf,
                var=var,
                fltr=fltr,
                dfr=dfr,
                sr=sr,
                stat.test=stat.test,
                group = group))
}

get.analysis <- function(IDs=NULL,group=NULL,save=FALSE,dataf=NULL){
    in.IT(c("plyr","dplyr","RCurl"))

    cat('# Load Key Table\n')
    # Load Key Table
    ML2.key <- get.GoogleSheet(data='ML2masteRkey')$df

    IDs=which(ML2.key$study.global.include==1)
    group  <- 'Source.Global'

    # Load Source Info Table
    ML2.sit <- get.GoogleSheet(url='https://docs.google.com/spreadsheets/d/1Qn_kVkVGwffBAmhAbpgrTjdxKLP1bb2chHjBMVyGl1s/export?format=csv')$df

    if(is.null(IDs)){IDs <- seq_along(ML2.key[,1])}

    # NOTE: CHANGE DIR TO FINAL LOCATION
    if(is.null(dataf)){
        cat('# Load data\n')
        ML2.S1 <- readRDS("~/Dropbox/Manylabs2/TestOutput/ML2_RawData_S1.rds")
        ML2.S2 <- readRDS("~/Dropbox/Manylabs2/TestOutput/ML2_RawData_S2.rds")
    }

    # Prepare list objects
    primary <- eval(parse(text= paste0('list(',paste0(unlist(ML2.key[IDs,'study.analysis']),collapse=' = list(), '),' = list())') ))
    data    <- eval(parse(text= paste0('list(',paste0(unlist(ML2.key[IDs,'study.analysis']),collapse=' = list(), '),' = list())') ))

    cat("# Start loop based on elements of IDs\n")
    for(s in IDs){
        cat("\n\n\n========",s,"===",ML2.key[[s,'study.analysis']],"========\n")
        cat("# Get the correct slate according to information in ML2.key['study.slate']\n")
        ifelse(ML2.key[s,'study.slate'] == 1, ML2.df <- ML2.S1,ML2.df<-ML2.S2)
        cat('# Add a unique ID as $uID\n')
        ML2.df$uID = seq(1,nrow(ML2.df))

        ML2 <- get.results(ML2.key,ML2.df,s,group)

        # Data list for calculating Effect Sizes CI based on NCP
        primary[[s]] <- list(analysis.id   = c(ML2.key[s,'study.id'],ML2.key[s,'study.slate'],ML2.key[s,'study.name']),
                             analysis.name = ML2.key[s,'study.analysis'],
                             analysis.group   = group,
                             stat.varfun   = ML2.key[s,'stat.vars'],
                             stat.type     = ML2.key[s,'stat.type'],
                             stat.ncp      = eval(parse(text=ML2.key[s,'stat.ncp'])),
                             stat.df       = eval(parse(text=ML2.key[s,'stat.df'])),
                             stat.N        = ML2$var$N,
                             stat.info     = ML2$info,
                             analysis.extra = list(stat.test = stat.test))

        # Data list for output to spreadsheet
        data[[s]] <- list(stat.analysis.name  = ML2.key[s,'study.analysis'],
                          stat.info           = ML2.in,
                          stat.data.cleanchain= ML2.id,
                          stat.data.raw       = ML2.sr$RawDataFilter,
                          stat.data.cleaned   = ML2.sr[1:length(ML2.sr)-1],
                          stat.data.analysed  = ML2.var,
                          stat.test.result    = stat.test)


        #     prime <- data.frame(cbind(study.id      = ML2.key[s,'study.id'],
        #                               study.slate   = ML2.key[s,'study.slate'],
        #                               study.name    = ML2.key[s,'study.name'],
        #                               analysis.name = ML2.key[s,'study.analysis'],
        #                               stat.varfun   = ML2.key[s,'stat.vars'],
        #                               stat.type     = ML2.key[s,'stat.type'],
        #                               stat.ncp      = NA, #eval(parse(text=ML2.key[s,'stat.ncp'])),
        #                               stat.df1      = NA, #eval(parse(text=ML2.key[s,'stat.df'])),
        #                               stat.df2      = NA,
        #                               stat.n1       = NA, #ML2$var$N,
        #                               stat.n2       = NA, #ML2$var$N,
        #                               analysis.group= NA)) #ML2$info)
        primer <- list()
        for(i in seq_along(ML2$var)){
            primed <- prime
            stat.test <- ML2$stat.test[[i]]
            #cat('# Data list for calculating Effect Sizes CI based on NCP\n')
            primed$stat.ncp <- eval(parse(text=ML2.key[s,'stat.ncp']))
            df <- eval(parse(text=ML2.key[s,'stat.df']))
            ifelse(length(df)<2,{df2 <- NA},{df2 <- df[2]})
            primed$stat.df1 <- df[1]
            primed$stat.df2 <- df2
            N <- ML2$var[[i]]$N[[1]]
            ifelse(length(N)<2,{n2 <- NA},{n2 <- N[2]})
            primed$stat.n1  <- N[1]
            primed$stat.n2  <- n2
            primed$analysis.group <- ML2$group[i]
            primer[[i]] <- primed
        }
        results <- ldply(primer)

        cat('# Data list for output to spreadsheet\n')
        data[[s]] <- list(stat.analysis.name  = ML2.key[s,'study.analysis'],
                          stat.info           = ML2$info,
                          stat.data.cleanchain= ML2$fltr,
                          stat.data.raw       = sapply(ML2$sr, function(rdat) rdat$RawDataFilter),
                          stat.data.cleaned   = sapply(ML2$sr, function(cdat) cdat[1:length(cdat)-1]),
                          stat.data.analysed  = ML2$var,
                          stat.test.result    = ML2$stat.test)
    }

    if(save){# Save to RData and xlsx
        setwd("~/Dropbox/Manylabs2/TestOutput")
        save.ML2(data,study=IDs)}

    #   x <- ML2.primary[IDs]
    #   y <- ML2.data[IDs]
    #   names(x) <- ML2.key[IDs,'study.analysis']
    #   names(y) <- ML2.key[IDs,'study.analysis']

    return(list(ML2.primary=primer,ML2.data=data))

}

# Save ML2 data -----------------------------------------------------------
save.ML2 <- function(data,type=c('all','csv','xlsx','R')[1],toDir=getwd(),prefix='ML2_data_',studies=1:length(data)){
    cat(toDir)
    switch(type,
           all  = toSave<-c(xlsx=TRUE ,R=TRUE ,csv=TRUE ),
           csv  = toSave<-c(xlsx=FALSE,R=FALSE,csv=TRUE ),
           xlsx = toSave<-c(xlsx=TRUE ,R=FALSE,csv=TRUE ),
           R    = toSave<-c(xlsx=FALSE,R=TRUE ,csv=FALSE)
    )
    outlist <- llply(seq_along(studies),function(s) (data[[s]]$stat.test.result))

    if(toSave['xlsx']){
        require(xlsx)
        # Save results for each analysis to spreadsheet, sources as
        for(s in studies){
            wb     <- createWorkbook(type="xlsx")
            rdf    <- melt(data=data[[s]]$stat.data.raw)

            sheet  <- createSheet(wb,sheetName=paste(unique(rdf$Source.Global)))
            addDataFrame(as.data.frame(rdf), sheet, startRow=1, startColumn=1)
            #autoSizeColumn(sheet,colIndex=1:ncol(data[[s]]$stat.data.precleaned[[tab]]))
            rm(sheet)

            #     # Save data for each  analysis to spreadsheet
            #     for(s in studies){
            #       if(length(sum(data[[s]]$stat.data.analysed$N))>0){
            #         wb     <- createWorkbook(type="xlsx")
            #
            #         sheet0 <-  createSheet(wb,sheetName="stat.test.result")
            #         addDataFrame(as.data.frame(capture.output(outlist[[s]])), sheet0, startRow=1, startColumn=1)
            #         rm(sheet0)
            #
            #         sheet1 <- createSheet(wb,sheetName="stat.data.analysed")
            #         header <- createRow(sheet1, 1)
            #         cols   <- createCell(header, colIndex=1:(length(data[[s]]$stat.data.analysed)-1))
            #         rows   <- createRow(sheet1, 2:(sum(data[[s]]$stat.data.analysed$N)+1))
            #         cells  <- createCell(rows, colIndex=1:(length(data[[s]]$stat.data.analysed)-1))
            #         for(c in 1:(length(data[[s]]$stat.data.analysed)-1)){
            #           if(is.data.frame(data[[s]]$stat.data.analysed[[c]])){
            #             addDataFrame(data[[s]]$stat.data.analysed[[c]], sheet1, startRow=1, startColumn=1,showNA=SN)
            #           } else {
            #             #addDataFrame(as.data.frame(data[[s]]$stat.data.analysed[[c]]), sheet1, startRow=1, startColumn=1,showNA=SN)
            #
            #             setCellValue(cols[[1,c]],names(data[[s]]$stat.data.analysed)[c],showNA=SN)
            #
            #             if(is.factor(data[[s]]$stat.data.analysed[[c]])){
            #               mapply(setCellValue, cells[seq(1,length(data[[s]]$stat.data.analysed[[c]])),c],paste(data[[s]]$stat.data.analysed[[c]]),showNA=SN)
            #             } else {
            #               mapply(setCellValue, cells[seq(1,length(data[[s]]$stat.data.analysed[[c]])),c],data[[s]]$stat.data.analysed[[c]],showNA=SN)
            #             }
            #           }
            #           #autoSizeColumn(sheet1,colIndex=c)
            #         }
            #         rm(sheet1,header,cols,rows,cells)
            #
            #         sheet2 <- createSheet(wb,sheetName="stat.data.raw")
            #         addDataFrame(data[[s]]$stat.data.raw, sheet2, startRow=1, startColumn=1,showNA=SN,characterNA="NA")
            #         #autoSizeColumn(sheet2,colIndex=1:ncol(data[[s]]$stat.data.raw))
            #         rm(sheet2)
            #
            #         #         for(tab in seq_along(data[[s]]$stat.data.cleaned)){
            #         #           sheet <- createSheet(wb,sheetName=paste("cleaned.",names(data[[s]]$stat.data.cleaned[tab])))
            #         #           addDataFrame(data[[s]]$stat.data.cleaned[[tab]], sheet, startRow=1, startColumn=1)
            #         #           #autoSizeColumn(sheet,colIndex=1:ncol(data[[s]]$stat.data.precleaned[[tab]]))
            #         #           rm(sheet)
            #         #         }
            #         saveWorkbook(wb,file=paste0(toDir,'/',prefix,data[[s]]$stat.analysis.name,".xlsx"))
        }
        saveWorkbook(wb,file=paste0(toDir,'/',prefix,data[[s]]$stat.analysis.name,Sys.Date(),".xlsx"))
    }

    else {
        cat('Skipped: ',s)
    }

    if(toSave['R']){
        save(data,file=paste0(toDir,'/',prefix,"ALL.RData"))
    }
}

save.ML2.results <- function(results,outlist,toDir=getwd(),prefix='ML2_results_',studies=1:nrow(results)){
    SN = FALSE
    require(xlsx)
    # Save results for each analysed to spreadsheet
    wb     <- createWorkbook(type="xlsx")
    sheet1 <- createSheet(wb,sheetName="ML2 Results Summary")
    addDataFrame(results, sheet1, startRow=1, startColumn=1,showNA=SN)
    for(s in studies){
        sheet <- createSheet(wb,sheetName=paste(names(outlist)[[s]]))
        addDataFrame(as.data.frame(capture.output(outlist[[s]])), sheet, startRow=1, startColumn=1)
        #autoSizeColumn(sheet,colIndex=1:ncol(data[[s]]$stat.data.precleaned[[tab]]))
        rm(sheet)
        saveWorkbook(wb,file=paste0(toDir,'/',prefix,Sys.Date(),".xlsx"))
    }
}

# varfuns - Functions used to prepare data for analysis ------------------------------

# func: Clean Source Label
clean.Source <- function(source.raw, SourceTable){
    source.clean <- gsub("([[:blank:]]|[[:punct:]]|[[:cntrl:]])+","",source.raw)
    for(s in seq_along(SourceTable$Source.Field.Raw)){
        ID <- which(source.clean%in%SourceTable$Source.Field.Raw[[s]])
        source.clean[ID] <- SourceTable$Source.Field.Correct[[s]]
        }
    return(as.data.frame(source.clean))
    }

# func: Get additional changes, variables, etc.
get.fieldAdd <- function(data,stable){
    data$Source.Global    <- NA
    data$Source.Primary   <- NA
    data$Source.Secondary <- NA
    data$Country      <- NA
    data$Language     <- NA
    data$Execution    <- NA
    data$SubjectPool  <- NA
    data$Setting      <- NA
    data$Tablet       <- NA
    data$Pencil       <- NA
    data$StudyOrder   <- NA
    data$IDiffOrder   <- NA

    for(s in seq_along(stable$Source)){
        ID <- which((stable$Source[[s]]==data$source)&(stable$Filename[[s]]==data$.id))
        if(length(ID)>0){
            data$Source.Global[ID]<- stable$Source.Global[[s]]
            data$Source.Primary[ID]   <- stable$Source.Global[[s]]
            data$Source.Secondary[ID] <- stable$Source[[s]]
            data$Country[ID]      <- stable$Country[[s]]
            data$Language[ID]     <- stable$Language[[s]]
            data$Execution[ID]    <- stable$Execution[[s]]
            data$SubjectPool[ID]  <- stable$SubjectPool[[s]]
            data$Setting[ID]      <- stable$Setting[[s]]
            data$Tablet[ID]       <- stable$Tablet[[s]]
            data$Pencil[ID]       <- stable$Pencil[[s]]
            data[ID, "StudyOrder"]   <- data[ID, stable$StudyOrder[[s]]]
            data[ID, "IDiffOrder"]   <- data[ID, stable$IDiffOrder[[s]]]
        }
    }
    return(as.data.frame(data))
}

scanColumns <- function(pattern, data){
    # Function to scan columns of a dataframe for a regex pattern and return a variable (in rows) by case (in columns) logical matrix
    # Includes some error catching, correcting and reporting
    idS    <- list()
    cNames <- colnames(data)
    for(c in seq(1,ncol(data))){
        tmp      <- try.CATCH(c(grepl(pattern,data[[c]],useBytes = T,ignore.case=T)))
        if(is.null(tmp$warning)){
            idS[[c]]<-tmp$value
        } else {
            cat('Error in column: ', colnames(data)[c],'->\t',paste0(tmp$warning),'\n')
        }
    }
    idD   <- ldply(idS)
    cs    <- rowSums(idD)
    colID <- which(cs>0)
    idS   <- t(idD[colID, ])

    return(list(idMatrix = idS,
                idColnames = cNames[colID])
    )
}

# func: Clean data fields to NA according to pattern
# Default: Set test trials and -99 to NA
clean.ML2fieldsNA <- function(source.raw, pattern="(test|-99)"){

    # First mark list of known test runs for removal
    disp(message = 'Marking known test session for removal', header = 'Clean ML2 Test Data - Step 1', footer = F)
    source.raw$Finished[which(c("R_blS3XHmD4p14kf3", "R_8Cjo1Lguez90bqt", "R_bC85NEijsCjOqep")%in%source.raw$ResponseID)] <- 0

    # Now search for 'test' and apply Finished == 1 filter
    disp(message = 'Checking columns "age", "sex", "hometown", "education" and "comments" for a variant of pattern: "test"', header = 'Clean ML2 Test Data - Step 2', footer = F)
    idS <- scanColumns(pattern='test',data = source.raw[, c('age','sex','hometown','education','comments')])
    if(sum(colSums(idS$idMatrix) > 0)){
        disp(message = paste0('Column "',idS$idColnames[colSums(idS$idMatrix)>0],'" contained cases with a variant of pattern: "test"'), header = F, footer = F)
        source.raw$Finished[(rowSums(idS$idMatrix)>0)] <- 0
    } else {
        source.clean <- source.raw
    }

    # Remove the test sessions
    source.clean <- source.raw %>% filter(Finished == 1)

    # Now look for '-99'
    disp(message = 'Checking all columns except "LocationLongitude" for pattern: "-99"', header = 'Clean Test ML2 Data - Step 3', footer = F)
    idS <- scanColumns(pattern='-99', data = source.clean[ ,which(!colnames(source.clean)%in%c("LocationLongitude"))])
    # Set -99 to NA
    for(c in seq(1,ncol(idS$idMatrix))){
        source.clean[idS$idMatrix[ ,c], as.numeric(colnames(idS$idMatrix)[c])] <- NA
    }
    disp(message = paste0('Column "',idS$idColnames[(colSums(idS$idMatrix)>0)],'" contained cases with pattern "-99"'), header = F, footer=F)
    disp(message = 'Done!')
    return(source.clean)
}

get.CSVdata <- function(files, finishedOnly=TRUE){
    if(grepl(".xlsx",files)){
        temp <- read.xlsx2(files, sheetIndex=1,stringsAsFactors=F)
        colnames(temp)[1:10] <- temp[1,1:10]
        df <- tbl_df(temp[-1,])
    } else {
        t1 <- scan(files,sep=",",nlines = 1,what="character", quiet = T)
        t2 <- scan(files,sep=",",skip = 1, nlines = 1,what="character", quiet = T)
        t1[1:10] <- t2[1:10]
        df <- try.CATCH(tbl_df(read.csv(files,skip=2,header=F,stringsAsFactors=F,col.names=t1)))
        ifelse(is.data.frame(df$value), {df <- df$value}, {df <- tbl_df(read.csv(files,skip = 2,header = F,stringsAsFactors = F,col.names = t1, fileEncoding="latin1"))})
    }
    if(finishedOnly){
        # Remove cases that did not complete the study
        df <- df %>% filter(Finished == 1)
    }
    return(df)
}

get.CSVinfo <- function(files){
    t1 <- scan(files,sep=",",nlines=1,what="character",quiet=T)
    t2 <- scan(files,sep=",",skip=1, nlines=1,what="character",quiet=T)
    df <- tbl_df(read.csv(files,skip=2,header=F,stringsAsFactors=F,col.names=t1))
    return(list(header1=t1,header2=t2,Info=data.frame(fileName=files,ncases=nrow(df),nvars=ncol(df))))
}

get.chain <- function(inf){
    # Build a filter chain

    filt.vars <- unlist(inf$study.cases.include,recursive=F)

    # include sites
    filt.site <- paste0(" %>% dplyr::filter(",paste0(inf$study.sites.include),")")
    #filt.site <- paste0(" %>% filter(is.character(.id))")

    # Data frame filter
    filt.df <- paste0(" %>% dplyr::select(", paste0(inf$id.vars,collapse=","),")", filt.site)

    #Variables filter
    return(list(df=filt.df,vars=filt.vars))
}

get.cases <- function(rule,study.vars,study.vars.labels,stat.params){
    #rule <- cases.include
    type  <-names(rule)
    if(!is.matrix(rule[[1]])){rule  <- rbind(rule[[1]])} else {rule <- rule[[1]]}
    Nrule <- nrow(rule)

    isna <- ifelse(stat.params$censorNA,{" & !is.na(X)"},{""})
    do <- list()

    # Rule 1 should always be about the variables in 'study.vars'
    # Rule 2..N can concern a subset of the variables in 'study.vars'. Assumption is: number of variables is a multiple of number of conditions, with variables grouped to fit to each condition.
    r = 1
    filtvars <- study.vars.labels[names(study.vars.labels)%in%rule[r,1]]

    # Start with 'pre'
    pre <-llply(unlist(filtvars), function(v){paste0(v,' %>% filter(')})

    if(all(filtvars[[1]]%in%names(study.vars))){
        filtvars <- study.vars
        do <- laply(seq_along(filtvars), function(v) laply(seq_along(filtvars[[v]]),
                                                           function(vv) paste0(gsub("X",filtvars[[v]][vv],rule[r,2]), gsub("X",filtvars[[v]][vv],isna))))
    }

    if(!is.matrix(do)){do <- as.matrix(do)}

    if(Nrule > 1){
        s <- ncol(do)
        for(r in 2:Nrule){
            s <- s+1
            filtvars <- study.vars.labels[names(study.vars.labels)%in%rule[r,1]]
            do  <- cbind(do, laply(seq_along(filtvars), function(vv) paste0(gsub("X",filtvars[[vv]],rule[r,2]))))
        }
    }
    case <- llply(seq_along(pre), function(fr) paste0(pre[[fr]], paste0(do[fr, ],collapse = ' & '),")"))
    names(case) <- names(study.vars)
    return(eval(parse(text=paste0('list(',rule[1],'=case)'))))
}

## OLD get.cases
# ifelse(rule[[1]][1]=="all",{fvars=unlist(study.vars)},{fvars=study.vars[names(study.vars)%in%unlist(study.vars.labels[rule[[1]][1]])]})
# Xrule.id <- grepl("X",rule[[1]])
#     # Fill in comands 'do'
#     case=list()
#     do <- llply(seq_along(fvars), function(v){
#         for(vi in seq(2,Nrule)){
#             if(Xrule.id[vi]){
#                 if(stat.params$censorNA){
#                     ifelse(type=="each",{
#                         case[[vi-1]] <- llply(fvars[[v]], function(vii){paste0(gsub("X",vii,rule[[1]][vi]),' & !is.na(',vii,')')})},{
#                             case[[vi-1]] <- paste0(llply(fvars[[v]], function(vii){paste0(gsub("X",vii,rule[[1]][vi]),' & !is.na(',vii,')')}), collapse=' & ')
#                         })
#                 } else {
#                     ifelse(type=="each",{
#                         case[[vi-1]] <- llply(fvars[[v]], function(vii){paste0(gsub("X",vii,rule[[1]][vi]))})},{
#                             case[[vi-1]] <- paste0(llply(fvars[[v]], function(vii){paste0(gsub("X",vii,rule[[1]][vi]))}), collapse=' & ')
#                         })
#                 }
#             } else {
#                 case[[vi-1]]<-paste0(rule[[1]][vi])
#             }
#         }
#         paste0(unlist(case),collapse=' & ')
#     })
#
#     case        <- llply(paste0(pre,do,')'))
#     names(case) <- names(fvars)

get.info <- function(keytable,cols){
    # Read Variables and Parameters from:
    #keytable <- ML2.key[s,]
    study.vars         <- eval(parse(text=keytable[,'study.vars']))
    study.vars.labels  <- eval(parse(text=keytable[,'study.vars.labels']))
    cases.include      <- eval(parse(text=keytable[,'study.cases.include']))
    stat.params        <- eval(parse(text=keytable[,'stat.params']))
    #cases <- llply(seq_along(cases.include),function(i) get.cases(cases.include[i],study.vars,study.vars.labels,stat.params))
    cases              <- get.cases(cases.include, study.vars, study.vars.labels, stat.params)
    sites.include      <- eval(parse(text=keytable[,'study.sites.include']))
    if(sites.include[[1]][1]=="all"){sites.include[[1]]<-'is.character(source)'}

    # Find correct columns in this dataset according to ML2.key: 'ML2.in$study.vars'
    id.vars  <- which(cols%in%c(unlist(study.vars),'uID','.id','age','sex','source','Source.Global','Source.Primary','Source.Secondary','Country','Language','SubjectPool','Setting','Tablet','Pencil','Execution', 'StudyOrderN','IDiffOrderN'))
    return(list(study.vars          = study.vars,
                study.vars.labels   = study.vars.labels,
                stat.params         = stat.params,
                study.cases.include = cases,
                study.sites.include = sites.include,
                id.vars             = id.vars))
}

get.sourceData <- function(ML2.id,ML2.df,ML2.in){
    N       <- numeric(length(ML2.in$study.vars))
    #study.vars[1]    <- unlist(ML2.in$study.vars)
    vars <- list()
    dfname <- list()
    RawData <- list()
    #id <- factor(ML2.df$.id)
    for(i in seq_along(ML2.in$study.vars)){
        dfname[i] <- names(ML2.in$study.vars)[[i]]
        eval(parse(text=paste0(names(ML2.in$study.vars)[[i]]," <- tbl_df(dplyr::select(ML2.df,", paste0(c(ML2.in$study.vars[[i]],'uID'),collapse=","),"))")))

        # Change text to numbers
        suppressWarnings(if(any(eval(parse(text= paste0('apply(',names(ML2.in$study.vars)[i],',2,is.numeric)==FALSE'))))){
            eval(parse(text=paste0(names(ML2.in$study.vars)[i],' <- data.frame(sapply(which(apply(',names(ML2.in$study.vars)[i],',2,is.numeric)==FALSE), function(si) as.numeric(',names(ML2.in$study.vars)[i],'[[si]])))')))
        })
    }

    for(i in seq_along(dfname)){
        eval(parse(text=paste0(dfname[i],' <- ', unlist(ML2.id$vars)[i])))
        N[i]         <- eval(parse(text=paste0("nrow(",dfname[i],")")))
        RawData[[i]] <- mutate(ML2.df, Included = eval(parse(text=paste0('ML2.df$uID %in% ',dfname[i],'$uID'))))
        eval(parse(text=paste(dfname[i],"<-", dfname[i]," %>% dplyr::select(which(colnames(",dfname[i],")!='uID'))")))
        eval(parse(text=paste0(dfname[i],' -> vars[[i]]')))
    }
    vars[[length(ML2.in$study.vars)+1]] <- N
    #if(length(ML2.in$study.vars.labels)==0){ML2.in$study.vars.labels <- list(NoLabels="None")}
    vars[[length(ML2.in$study.vars)+2]] <- ML2.in$study.vars.labels
    vars[[length(ML2.in$study.vars)+3]] <- RawData
    names(vars) <- c(names(ML2.in$study.vars),"N","labels","RawDataFilter")
    return(vars)
}

get.Stats <- function(indata){
    # Search patterns
    real_pat <- paste("(\\d*(","(\\,|\\.)?","\\d+){0,3})",sep="[[:blank:]]?")

    type.Chi2  <- paste("(?<type>\\b(([xXχ])|([cC]hi)|([\\xcf\\x87]))+","[.^]?","(2|(sq[uare]))*)",sep="[[:blank:]]*")
    type.t     <- paste("(?<type>\\b[tT]+",")",sep="[[:blank:]]*")
    type.F     <- paste("(?<type>\\b[fF]+(1|2)?",")",sep="[[:blank:]]*")
    type.r     <- paste("(?<type>\\b[rRρΡ]+",")",sep="[[:blank:]]*")
    type.Z     <- paste("(?<type>\\b[zZ]+",")",sep="[[:blank:]]*")
    type.B     <- paste("(?<type>\\b[bB]+","[s]?",")",sep="[[:blank:]]*")
    type.beta  <- paste("(?<type>\\b([ββ]|[bB]eta)+","[hat]?",")",sep="[[:blank:]]*")
    type.pr    <- paste("(?<type>\\b(pr))+",sep="[[:blank:]]*")

    statsPatterns <- list(
        F    = paste(type.F,"[(](?<df1>",real_pat,")[,]","(?<df2>",real_pat,")","[)]","=","(?<stat>",real_pat,")",sep="[[:blank:]]*"),
        t    = paste(type.t,"[(](?<df1>",real_pat,")[)]","=","(?<stat>","[-]?",real_pat,")",sep="[[:blank:]]*"),
        X2   = paste(type.Chi2,"[(](?<df1>",real_pat,")","([,]","N","=","(?<N>",real_pat,"))?[)]","=(?<stat>",real_pat,")",sep="[[:blank:]]*"),
        r    = paste(type.r,"([(](?<df1>",real_pat,")[)])?","=(?<stat>","[-]?",real_pat,")",sep="[[:blank:]]*"),
        Z    = paste(type.Z,"=(?<stat>",real_pat,")",sep="[[:blank:]]*"),
        B    = paste(type.B,"=(?<stat>",real_pat,")",sep="[[:blank:]]*"),
        beta = paste(type.beta,"=(?<stat>",real_pat,")",sep="[[:blank:]]*"),
        pr   = paste(type.pr,"=(?<stat>",real_pat,")",sep="[[:blank:]]*")
    )
    #data<-iconv(data,from="CP1252",to="UTF-8",mark=F)
    data    <- gsub("[[:blank:]]+","", indata, perl=T)
    outdata <- data.frame(rows=1:length(data),cbind(stat.type=NA,stat.ncp=NA,stat.df1=NA,stat.df2=NA))
    IDs.ori <- llply(statsPatterns,function(id) grep(id, data, ignore.case = T, perl=T))
    cleanPat <- "([[:blank:]]|=)*"
    for(m in seq_along(data)){
        patID <- sapply(IDs.ori, function(idp) m%in%idp)
        if(any(patID)){
            matchOri <- gregexpr(pattern=statsPatterns[patID],text=data[m],ignore.case=TRUE,perl=TRUE)[[1]]
            for(stat.field in attributes(matchOri)$capture.names[attributes(matchOri)$capture.names!='']){
                fieldID <- attributes(matchOri)$capture.names==stat.field
                if(attributes(matchOri)$capture.start[fieldID]!=-1){
                    switch(stat.field,
                           type = outdata$stat.type[m] <- names(statsPatterns[patID]),
                           stat = outdata$stat.ncp[m] <- gsub(cleanPat,"",(substr(data[[m]], start=attributes(matchOri)$capture.start[fieldID],stop=attributes(matchOri)$capture.start[fieldID]+attributes(matchOri)$capture.length[fieldID]-1))),
                           df1 = outdata$stat.df1[m] <- gsub(cleanPat,"",(substr(data[[m]], start=attributes(matchOri)$capture.start[fieldID],stop=attributes(matchOri)$capture.start[fieldID]+attributes(matchOri)$capture.length[fieldID]-1))),
                           df2 = outdata$stat.df2[m] <- gsub(cleanPat,"",(substr(data[[m]], start=attributes(matchOri)$capture.start[fieldID],stop=attributes(matchOri)$capture.start[fieldID]+attributes(matchOri)$capture.length[fieldID]-1)))
                           #N = outdata$stat.N[m] <- gsub(cleanPat,"",(substr(data[[m]], start=attributes(matchOri)$capture.start[fieldID],stop=attributes(matchOri)$capture.start[fieldID]+attributes(matchOri)$capture.length[fieldID]-1)))
                    )
                }
            }
        } else {
            outdata$stat.type[m]=paste("unknown:",regmatches(data[m],matchOri))
        }
    }
    return(outdata)
}

# STUDY variable functions [ML2] ------------------------------------------------

# Huang.1
varfun.Huang.1 <- function(vars=ML2.sr){
    # Analysis plan. The coordinates of the click on the map will be recorded (X, Y) from the
    # top-left of the image. The mean difference between the high and low-SES conditions for
    # north/south location of click (Y) will be compared with an independent samples t-test. All
    # participants who indicate an area within the boundaries of the map will be included in the
    # analysis.
    # [...]

    # huan1.1_Y1 = Y position of the mouse (high ses condition).
    # huan2.1_Y1 = Y position of the mouse (low ses).
    # smaller numbers=upper position.

    # huan1.1_R0 and huan2.1_R0 indicate for each condition whether they clicked inside the map (1) or outside (0).
    # For each condition they must have clicked inside the map (=1) to be included in the analysis.


    # ** Huang.1) Whole group comparing SESvignette hi vs lo (global and primary)
    # Huang.2) Whole group disregarding SESvignette comparing Homewealth north vs south (secondary)
    # Huang.3) only hongkong disregarding Homewealth comparing SESvignette hi vs lo (secondary)
    # Huang.4) only US disregarding Homewealth comparing SESvignette hi vs lo (secondary)
    # Huang.5) only Tablets & homewealth = north comparing SESvignette hi vs lo (secondary)

    return(list(High = vars$High[[1]],
                Low  = vars$Low[[1]],
                N    = vars$N))
}

varfun.Huang.2 <- function(vars=ML2.sr){
    # Analysis plan.
    #
    # [...]
    #
    # The test for replicating the cultural difference observed in Huang et al. will be conducted
    # on a subset of the participants that respond on the wealth in hometown individual difference item
    # that wealthy people tend to live in the North (akin to original U.S. sample) versus wealthy people
    # tend to live in the South (akin to original Hong Kong sample). The entire sample will be used
    # for investigating variation in effects across sample and setting.

    # huan1.1_Y1 = Y position of the mouse (high ses condition).
    # huan2.1_Y1 = Y position of the mouse (low ses).
    # smaller numbers =upper position

    # Huang.1) Whole group comparing SESvignette hi vs lo (global and primary)
    # ** Huang.2) Whole group disregarding SESvignette comparing Homewealth north vs south (secondary)
    # Huang.3) only hongkong disregarding Homewealth comparing SESvignette hi vs lo (secondary)
    # Huang.4) only US disregarding Homewealth comparing SESvignette hi vs lo (secondary)
    # Huang.5) only Tablets & homewealth = north comparing SESvignette hi vs lo (secondary)

    return(list(North = c(vars$High[[1]],vars$Low[[1]])[factor(c(vars$High[[3]],vars$Low[[3]]),levels=c(1,2),vars$labels$Response)=='North'],
                South = c(vars$High[[1]],vars$Low[[1]])[factor(c(vars$High[[3]],vars$Low[[3]]),levels=c(1,2),vars$labels$Response)=='South'],
                N    = vars$N))
}

varfun.Kay.1 <- function(vars=ML2.sr){
    # Analysis plan. Following Kay et al. (2014), we will create an index of willingness to engage in goal pursuit for each participant by
    # (1) regressing the mean of the two goal pursuit items on the centered mean of the goal subjective value item,
    # (2) calculating the unstandardized residual for each participant, and
    # (3) add to those the mean value for the self-regulation items measuring willingness to engage in goal pursuit.
    # Then, the two conditions will be compared using an independent samples t-test.

    var.Order <- rowMeans(select(vars$Order,kay1.5,kay1.6))
    # [(1) Centered Subjective value = (GP1 + GP2)]$(2)residuals + (3) mean Willingness to engage in goal pursuit
    var.windex.Order <-lm(scale(select(vars$Order,kay1.4),scale=F)~var.Order)$residuals + colMeans(select(vars$Order,kay1.5))

    var.DisOrder <- rowMeans(select(vars$DisOrder,kay2.5,kay2.6))
    # [(1) Centered Subjective value = (GP1 + GP2)]$(2)residuals + (3) mean Willingness to engage in goal pursuit
    var.windex.DisOrder <-lm(scale(select(vars$DisOrder,kay2.4),scale=F)~var.DisOrder)$residuals + colMeans(select(vars$DisOrder,kay2.5))

    return(list(Order=var.windex.Order,
                DisOrder=var.windex.DisOrder,
                N=vars$N))
}

varfun.Alter.1 <- function(vars=ML2.sr){
    # Analysis plan. Similar to Alter et al. (2007), we will conduct an independent samples ttest to determine whether accuracy in solving moderately difficult syllogisms differ by font condition (fluent versus disfluent). The original study focused on the moderately difficult questions, on the basis that participants’ performance could vary enough to detect changes in processing depth. Our primary analysis strategy will be sensitive to potential differences across samples in ability on syllogisms. We will first determine which syllogisms were moderately difficult to participants by excluding any of the six items, within each sample, that were answered correctly by fewer than 25% of participants or more than 75% of participants across conditions. The remaining syllogisms will be the basis of computing mean syllogism performance for each participant. For a direct comparison with the original effect size, only English in-lab samples will be used for two reasons: (1) we cannot adequately control for online participants “zooming in” on the page or otherwise making the font more readable, and (2) a different font may be used in some translated versions because the original font (Myriad Web) may not support all languages. All samples will be included in the investigation of cross-site variability in effect size.

    var.correct <- list(s1=c(7),
                        s2=c(8),
                        s3=c(5,6,7,8),
                        s4=c(8),
                        s5=c(3),
                        s6=c(8))

    lowP <- .25
    hiP  <- .75

    # Get correct answers
    ok.Fluent   <- sapply(seq_along(vars$Fluent), function(c) unlist(vars$Fluent[,c])%in%var.correct[[c]])
    ok.DisFluent<- sapply(seq_along(vars$DisFluent), function(c) unlist(vars$DisFluent[,c])%in%var.correct[[c]])

    # Find columns
    # Syllogisms to include for each sample
    # INCLUSION PERCENTAGE BASED ON
    # FLUENT / DISFLUENT SEPERATELY: 1 5 6
    # BOTH: 1 5 6
    id.Fluent.cols    <- which((colSums(ok.Fluent)/nrow(ok.Fluent)>lowP)&(colSums(ok.Fluent)/nrow(ok.Fluent)<hiP))
    id.DisFluent.cols <- which((colSums(ok.DisFluent)/nrow(ok.DisFluent)>lowP)&(colSums(ok.DisFluent)/nrow(ok.DisFluent)<hiP))

    return(list(Fluent   = laply(1:length(cbind(ok.Fluent[,id.Fluent.cols])[,1]), function(c) sum(ok.Fluent[c,id.Fluent.cols])),
                DisFluent= laply(1:length(cbind(ok.DisFluent[,id.DisFluent.cols])[,1]), function(c) sum(ok.DisFluent[c,id.Fluent.cols])),
                N=vars$N,
                SylCorrect=list(Fluent=colSums(ok.Fluent)/nrow(ok.Fluent),
                                DisFluent=colSums(ok.DisFluent)/nrow(ok.DisFluent))
    )
    )
}

varfun.Alter.2 <- function(vars=ML2.sr){
    # Analysis plan. Similar to Alter et al. (2007), we will conduct an independent samples ttest to determine whether accuracy in solving moderately difficult syllogisms differ by font condition (fluent versus disfluent). The original study focused on the moderately difficult questions, on the basis that participants’ performance could vary enough to detect changes in processing depth. Our primary analysis strategy will be sensitive to potential differences across samples in ability on syllogisms.

    # As a secondary analysis, we will use the same two syllogisms from Alter et al (2007) for analysis regardless of performance to perfectly mirror the original analysis.

    var.correct <- list(s1=c(7),
                        s2=c(8),
                        s3=c(5,6,7,8),
                        s4=c(8),
                        s5=c(3),
                        s6=c(8))

    # Get correct answers
    ok.Fluent   <- sapply(seq_along(vars$Fluent), function(c) unlist(vars$Fluent[,c])%in%var.correct[[c]])
    ok.DisFluent<- sapply(seq_along(vars$DisFluent), function(c) unlist(vars$DisFluent[,c])%in%var.correct[[c]])

    # Syllogisms to include for each sample
    # First and last
    id.Fluent.cols      <- c(1,6)
    id.DisFluent.cols   <- c(1,6)

    return(list(Fluent    = rowSums(ok.Fluent[,id.Fluent.cols]),
                DisFluent = rowSums(ok.DisFluent[,id.DisFluent.cols]),
                N=vars$N))
}


varfun.Alter.3 <- function(vars=ML2.sr){
    # Analysis plan. Similar to Alter et al. (2007), we will conduct an independent samples ttest to determine whether accuracy in solving moderately difficult syllogisms differ by font condition (fluent versus disfluent). The original study focused on the moderately difficult questions, on the basis that participants’ performance could vary enough to detect changes in processing depth. Our primary analysis strategy will be sensitive to potential differences across samples in ability on syllogisms.

    # As a secondary analysis, we will use the same two syllogisms from Alter et al (2007) for analysis regardless of performance to perfectly mirror the original analysis.
    #
    # This function selects only cases for which Alter was the first study.

    var.correct <- list(s1=c(7),
                        s2=c(8),
                        s3=c(5,6,7,8),
                        s4=c(8),
                        s5=c(3),
                        s6=c(8))

    # Get ids for Alter first
    id <- sapply(seq_along(vars$RawDataFilter[[1]]$.id), function(i) unlist(strsplit(x = vars$RawDataFilter[[1]]$StudyOrderN[i], split = "[|]"))[[1]] == "Alter")

     # Get correct answers for ids Alter first
    if(sum(id,na.rm=T)>0){
        ok.Fluent   <- rbind(sapply(seq_along(vars$Fluent), function(c) unlist(vars$Fluent[id, c])%in%var.correct[[c]]))
        ok.DisFluent<- rbind(sapply(seq_along(vars$DisFluent), function(c) unlist(vars$DisFluent[id, c])%in%var.correct[[c]]))
   } else {
        ok.Fluent    <- rep(FALSE,6)
        ok.DisFluent <- rep(FALSE,6)
    }

    # Use 1,5,6
    id.Fluent.cols    <- c(1,5,6)
    id.DisFluent.cols <- c(1,5,6)

    return(list(Fluent    = rowSums(rbind(ok.Fluent[ ,id.Fluent.cols])),
                DisFluent = rowSums(rbind(ok.DisFluent[ ,id.DisFluent.cols])),
                N = vars$N,
                SylCorrect = list(Fluent=colSums(ok.Fluent)/nrow(ok.Fluent),
                                  DisFluent=colSums(ok.DisFluent)/nrow(ok.DisFluent)))
           )
}


varfun.Alter.4 <- function(vars=ML2.sr){
    # Analysis plan. Similar to Alter et al. (2007), we will conduct an independent samples ttest to determine whether accuracy in solving moderately difficult syllogisms differ by font condition (fluent versus disfluent). The original study focused on the moderately difficult questions, on the basis that participants’ performance could vary enough to detect changes in processing depth. Our primary analysis strategy will be sensitive to potential differences across samples in ability on syllogisms.

    # As a secondary analysis, we will use the same two syllogisms from Alter et al (2007) for analysis regardless of performance to perfectly mirror the original analysis.

    var.correct <- list(s1=c(7),
                        s2=c(8),
                        s3=c(5,6,7,8),
                        s4=c(8),
                        s5=c(3),
                        s6=c(8))

    # Get ids for Alter first
    id <- sapply(seq_along(vars$RawDataFilter[[1]]$.id), function(i) unlist(strsplit(x = vars$RawDataFilter[[1]]$StudyOrderN[i], split = "[|]"))[[1]] == "Alter")

  # Get correct answers for ids Alter first
    if(sum(id,na.rm=T)>0){
        ok.Fluent   <- rbind(sapply(seq_along(vars$Fluent), function(c) unlist(vars$Fluent[id, c])%in%var.correct[[c]]))
        ok.DisFluent<- rbind(sapply(seq_along(vars$DisFluent), function(c) unlist(vars$DisFluent[id, c])%in%var.correct[[c]]))
   } else {
        ok.Fluent    <- rep(FALSE,6)
        ok.DisFluent <- rep(FALSE,6)
    }

    # Syllogisms to include for each sample
    # First and last
    id.Fluent.cols      <- c(1,6)
    id.DisFluent.cols   <- c(1,6)

    return(list(Fluent    = rowSums(rbind(ok.Fluent[ ,id.Fluent.cols])),
                DisFluent = rowSums(rbind(ok.DisFluent[ ,id.DisFluent.cols])),
                N = vars$N,
                SylCorrect = list(Fluent=colSums(ok.Fluent)/nrow(ok.Fluent),
                                  DisFluent=colSums(ok.DisFluent)/nrow(ok.DisFluent)))
           )
}


varfun.Graham.1 <- function(vars=ML2.sr){
    # The 6 items for harm and fairness will be averaged to create a
    # “individualizing” foundations moral relevance score and the 9 items for ingroup, authority, and
    # purity will be averaged to create a “binding” foundations moral relevance score. The
    # relationship between political ideology and the “binding” and “individualizing” aggregates will
    # be calculated using zero order correlations.

    # The primary target of replication is the relationship
    # of political ideology with the “binding” foundations, and the relationship of political ideology
    # with the “individualizing” foundations is a secondary replication. All participants who complete
    # the corresponding measures will be included in analysis.

    return(list(Politics = vars$Binding$politics,
                Binding  = rowMeans(select(vars$Binding,-politics),na.rm=T),
                N        = vars$N))
}

varfun.Graham.2 <- function(vars=ML2.sr){
    # The 6 items for harm and fairness will be averaged to create a
    # “individualizing” foundations moral relevance score and the 9 items for ingroup, authority, and
    # purity will be averaged to create a “binding” foundations moral relevance score. The
    # relationship between political ideology and the “binding” and “individualizing” aggregates will
    # be calculated using zero order correlations. The primary target of replication is the relationship
    # of political ideology with the “binding” foundations, and the relationship of political ideology
    # with the “individualizing” foundations is a secondary replication. All participants who complete
    # the corresponding measures will be included in analysis.

    return(list(Politics   = vars$Individual$politics,
                Individual = rowMeans(select(vars$Individual,-politics)),
                N          = vars$N))
}

varfun.Rottenstreich.1 <- function(vars=ML2.sr){
    # Analysis plan. A two-way contingency table will be built with certainty condition (lowprobability
    # vs. certain) and choice (monetary reward vs. meeting favorite movie star) as factors.
    # The critical replication hypothesis will be given by a χ2 test and the effect size by an odds ratio.
    # All participants with valid data on the response will be included in the analysis.
    return(list(Choice = factor(c(vars$Certain[[1]],vars$Low[[1]]),levels=c(1,2),labels=vars$labels$Condition),
                Chance = factor(c(rep(1,vars$N[1]),rep(2,vars$N[2])),levels=c(1,2),labels=vars$labels$Response),
                N      = vars$N))
}

varfun.Bauer.1 <- function(vars=ML2.sr){
    # Analysis plan. We will compare the mean trust levels between conditions with an
    # independent samples t-test. All participants with data will be included in analysis
    #
    # Known differences from original. The original experiment included four additional
    # dependent variables: (1) responsibility for the crisis, (2) obligation to cut water usage, (3) how
    # much they viewed others as partners, and (4) how much others should use less water. The central
    # replication will be on the trust variable, while the other four dependent variables will be retained
    # in the procedure but not analyzed for the focal replication.

    return(list(Consumer  = vars$Consumer[[2]],
                Individual= vars$Individual[[2]],
                N         = vars$N))
}

varfun.Miyamoto.1 <- function(vars=ML2.sr){
    # Analysis plan. An ANCOVA will compare the mean estimates of the author’s true attitude across the two conditions, covarying for perceived constraint.
    # miya1.5=true attitude (pro-death condition; higher values=higher support for death penalty);
    # miya1.7=perceived constraint (pro-death condition; higher values=higher freedom);
    #
    # miya2.5=true attitude (against death penalty condition; higher values=higher support for death penalty);
    # miya2.7=perceived constraint (against death condition; higher values= higher freedom).

    return(list(Attitude  = c(vars$CapitalPro[[1]],vars$CapitalCon[[1]]),
                Constraint= c(vars$CapitalPro[[2]],vars$CapitalCon[[2]]),
                Condition = factor(c(rep(1,vars$N[1]),rep(2,vars$N[2])),levels=c(1,2),labels=vars$labels$Condition),
                N = vars$N))
}

varfun.Miyamoto.2 <- function(vars=ML2.sr){
    # Analysis plan. An ANCOVA will compare the mean estimates of the author’s true attitude across the two conditions, covarying for perceived constraint.
    # miya1.5=true attitude (pro-death condition; higher values=higher support for death penalty);
    # miya1.7=perceived constraint (pro-death condition; higher values=higher freedom);
    #
    # miya2.4=true attitude (against death penalty condition; higher values=higher support for death penalty);
    # miya2.7=perceived constraint (against death condition; higher values= higher freedom).

    return(list(Attitude  = c(vars$CapitalPro[[1]],vars$CapitalCon[[1]]),
                Constraint= c(vars$CapitalPro[[2]],vars$CapitalCon[[2]]),
                Condition = factor(c(rep(1,vars$N[1]),rep(2,vars$N[2])),levels=c(1,2),labels=vars$labels$Condition),
                N = vars$N))
}


varfun.Inbar.1 <- function(vars=ML2.sr){
    # Analysis plan. The five items of the contamination subscale of the Disgust Sensitivity
    # Scale-Revised will be averaged to create a single index of disgust sensitivity. For the primary
    # analysis, we will compute a correlation between disgust sensitivity and assessments of the
    # director’s intentionality in both the gay kissing and kissing conditions, and then compare the
    # correlations with an r-to-z transformation. The other two outcome measures will be examined as
    # secondary analyses following the same analysis strategy. All participants with relevant data will
    # be included in analysis.

    # disg1.11,disg1.12,disg2.10,disg2.12,disg2.13; --responses on the DS-R are scored as follows: True 1, False 0; Not disgusting 0, Slightly disgusting 0.5, Very disgusting 1

    vars$SameKiss$disg1.11 <- -vars$SameKiss$disg1.11
    vars$SameKiss$disg1.12 <- -vars$SameKiss$disg1.12
    vars$DiffKiss$disg1.11 <- -vars$DiffKiss$disg1.11
    vars$DiffKiss$disg1.12 <- -vars$DiffKiss$disg1.12

    vars$SameKiss <- mutate(vars$SameKiss,DSRs=rowMeans(scale.R(select(vars$SameKiss, starts_with("disg")))))
    vars$DiffKiss <- mutate(vars$DiffKiss,DSRd=rowMeans(scale.R(select(vars$DiffKiss, starts_with("disg")))))

    outcome <- c("Intent","Wrong","Encourage")[colnames(vars$SameKiss)[1]==c("inba1.3","inba1.4","inba1.5")]

    colnames(vars$SameKiss)[1]     <- outcome
    colnames(vars$DiffKiss)[1] <- outcome

    return(list(r1=cbind(vars$SameKiss['DSRs'],vars$SameKiss[outcome]),
                r2=cbind(vars$DiffKiss['DSRd'],vars$DiffKiss[outcome]),
                N = vars$N)
           )
}


varfun.Inbar.2 <- function(vars=ML2.sr){
    # Analysis plan. The five items of the contamination subscale of the Disgust Sensitivity
    # Scale-Revised will be averaged to create a single index of disgust sensitivity. For the primary
    # analysis, we will compute a correlation between disgust sensitivity and assessments of the
    # director’s intentionality in both the gay kissing and kissing conditions, and then compare the
    # correlations with an r-to-z transformation. The other two outcome measures will be examined as
    # secondary analyses following the same analysis strategy. All participants with relevant data will
    # be included in analysis.

    # disg1.11,disg1.12,disg2.10,disg2.12,disg2.13; --responses on the DS-R are scored as follows: True 1, False 0; Not disgusting 0, Slightly disgusting 0.5, Very disgusting 1

    vars$SameKiss$disg1.11 <- -vars$SameKiss$disg1.11
    vars$SameKiss$disg1.12 <- -vars$SameKiss$disg1.12
    vars$DiffKiss$disg1.11 <- -vars$DiffKiss$disg1.11
    vars$DiffKiss$disg1.12 <- -vars$DiffKiss$disg1.12

    vars$SameKiss <- mutate(vars$SameKiss,DSRs=rowMeans(scale.R(select(vars$SameKiss, starts_with("disg")))))
    vars$DiffKiss <- mutate(vars$DiffKiss,DSRd=rowMeans(scale.R(select(vars$DiffKiss, starts_with("disg")))))

    return(list(Fisher = list(r1=cor(vars$SameKiss['DSRs'],vars$SameKiss[[1]],use="pairwise.complete.obs"),
                              r2=cor(vars$DiffKiss['DSRd'],vars$DiffKiss[[1]],use="pairwise.complete.obs"),
                              n1=vars$N[1],
                              n2=vars$N[2]),
                N = vars$N)
    )
}


varfun.Inbar.3 <- function(vars=ML2.sr){
    # Analysis plan. The five items of the contamination subscale of the Disgust Sensitivity
    # Scale-Revised will be averaged to create a single index of disgust sensitivity. For the primary
    # analysis, we will compute a correlation between disgust sensitivity and assessments of the
    # director’s intentionality in both the gay kissing and kissing conditions, and then compare the
    # correlations with an r-to-z transformation. The other two outcome measures will be examined as
    # secondary analyses following the same analysis strategy. All participants with relevant data will
    # be included in analysis.

    # disg1.11,disg1.12,disg2.10,disg2.12,disg2.13; --responses on the DS-R are scored as follows: True 1, False 0; Not disgusting 0, Slightly disgusting 0.5, Very disgusting 1

    vars$SameKiss$disg1.11 <- -vars$SameKiss$disg1.11
    vars$SameKiss$disg1.12 <- -vars$SameKiss$disg1.12
    vars$DiffKiss$disg1.11 <- -vars$DiffKiss$disg1.11
    vars$DiffKiss$disg1.12 <- -vars$DiffKiss$disg1.12

    vars$SameKiss <- mutate(vars$SameKiss,DSRs=rowMeans(scale.R(select(vars$SameKiss, starts_with("disg")))))
    vars$DiffKiss <- mutate(vars$DiffKiss,DSRd=rowMeans(scale.R(select(vars$DiffKiss, starts_with("disg")))))

    return(list(Fisher = list(r1=cor(vars$SameKiss['DSRs'],vars$SameKiss[[1]],use="pairwise.complete.obs"),
                              r2=cor(vars$DiffKiss['DSRd'],vars$DiffKiss[[1]],use="pairwise.complete.obs"),
                              n1=vars$N[1],
                              n2=vars$N[2]),
                N = vars$N)
    )
}

varfun.Critcher.1 <- function(vars=ML2.sr){
    # crit1.1= % sold P97 in the US;
    # crit2.1= % sold P17 in the US
    # Analysis plan. The means for the P97 and P17 groups will be compared with an
    # independent samples t-test. Participants whose answers will not fall between 0 and 100 will be excluded from analysis.
    #   df.P97 <- select(tbl_df(ML2.df),which(colnames(ML2.df)%in%ML2.in$study.vars$Condition[1]))
    #   df.P97 <- slice(df.P97, which((ML2.id[[1]][,1]==T)&(ML2.id[[2]][,1]==T)))
    #
    #   df.P17 <- select(tbl_df(ML2.df),which(colnames(ML2.df)%in%ML2.in$study.vars$Condition[2]))
    #   df.P17 <- slice(df.P17, which((ML2.id[[1]][,2]==T)&(ML2.id[[2]][,2]==T)))
    #
    #   id.P97 <- ML2.df[ML2.id[[1]][,1]&ML2.id[[2]][,1],ML2.in$study.vars$Condition[1]]
    #   id.P17 <- ML2.df[ML2.id[[1]][,2]&ML2.id[[2]][,2],ML2.in$study.vars$Condition[2]]
    #
    return(list(P97    = as.numeric(vars$P97[[1]]),
                P17    = as.numeric(vars$P17[[1]]),
                N      = vars$N))
}

varfun.vanLange.1 <- function(vars=ML2.sr){
    # van.p.1.2_1 TO van.p.1.2_6 are the items of SVO measure.
    #
    # murphy et al. (2011) scoring: SVO degress=arctan [(mean Alloc other − 50)/(mean Allocation self - 50)].
    #
    # See SVO codes (this doc) for the list of paired amounts.
    # van.p2.1_1_TEXT= # of older siblings;
    # van.p2.1_2_TEXT=# of younger siblings.
    #
    # Analysis plan. The current replication focuses only on the observed direct positive
    # correlation between greater prosocial orientation and number of siblings. Participants must
    # respond to all 6 items of the SVO slider and have valid data for the two questions asking about
    # older and younger siblings, to be included in the analysis. Total number of siblings will be
    # obtained by adding the number of younger and older siblings. SVO slider scores will be scored
    # in according to the procedure recommended by Murphy et al. (2011). The resulting SVO slider
    # score will be correlated with the total number of siblings for the critical test.

    SVO <-rbind(cbind(c( 85,85 ), c( 85,76 ), c( 85,68 ), c( 85,59 ), c( 85,50 ), c( 85,41 ), c( 85,33 ), c( 85,24 ), c( 85,15 )),
                cbind(c( 85,15 ), c( 87,19 ), c( 89,24 ), c( 91,28 ), c( 93,33 ), c( 94,37 ), c( 96,41 ), c( 98,46 ), c( 100,50 )),
                cbind(c( 50,100 ), c( 54,98 ), c( 59,96 ), c( 63,94 ), c( 68,93 ), c( 72,91 ), c( 76,89 ), c( 81,87 ), c( 85,85 )),
                cbind(c( 50,100 ), c( 54,89 ), c( 59,79 ), c( 63,68 ), c( 68,58 ), c( 72,47 ), c( 76,36 ), c( 81,26 ), c( 85,15 )),
                cbind(c( 100,50 ), c( 94,56 ), c( 88,63 ), c( 81,69 ), c( 75,75 ), c( 69,81 ), c( 63,88 ), c( 56,94 ), c( 50,100 )),
                cbind(c( 100,50 ), c( 98,54 ), c( 96,59 ), c( 94,63 ), c( 93,68 ), c( 91,72 ), c( 89,76 ), c( 87,81 ), c( 85,85 )))


    SVO.self  <- SVO[seq(1,11,by=2),]
    SVO.other <- SVO[seq(2,12,by=2),]


    SVO.index <- ldply(seq(1,vars$N[1]), function(s){atan(
        (mean(SVO.other[array(c(1:6,unlist(vars$SVO[s, ])),dim=c(6,2))])-50)/
            (mean( SVO.self[array(c(1:6,unlist(vars$SVO[s, ])),dim=c(6,2))])-50))
    })

    vars$SVO$Siblings <- vars$SVO$van.p2.1_1_TEXT+vars$SVO$van.p2.1_2_TEXT

    return(list(SVO.index = SVO.index$V1,
                Younger   = vars$SVO$van.p2.1_1_TEXT,
                Older     = vars$SVO$van.p2.1_2_TEXT,
                Siblings  = vars$SVO$Siblings,
                N         = vars$N))
}

varfun.Hauser.1 <- function(vars=ML2.sr){
    #   haus1.1t = timing (side effect scenario);
    #   haus2.1t = timing (greater good scenario);
    #   haus1.2=previous experience (drop if 1 (yes));
    #   haus2.2=previous experience (drop if 1 (yes));
    #   haus1.1=morally permissible (side effect scenario; Yes=1);
    #   haus2.1=morally permissible (greater good scenario; Yes=1).

    # Analysis plan. Subjects will be excluded from all analyses if they take fewer than four
    # seconds to read and respond to either of the target scenarios. For the key confirmatory test
    # comparing with the original effect size, we will include only participants that indicate having no
    # prior experience with the task. The original authors suggested the effect may be weaker in
    # participants with prior exposure. Prior exposure will be investigated as a moderator for the other
    # analyses. A two-way contingency table will be built with Scenario (Means vs. Side-effect) and
    # Response (Yes vs. No) as factors. The critical replication hypothesis will be given by a one tailed
    # chi square test and the effect size by an odds ratio.
#     SideEffect  <- filter(vars$SideEffect, vars$SideEffect[vars$labels$Experience[1]]!=1 & vars$SideEffect[vars$labels$Timing[1]]>4)
#     GreaterGood <- filter(vars$GreaterGood,vars$GreaterGood[vars$labels$Experience[2]]!=1 & vars$GreaterGood[vars$labels$Timing[2]]>4)
     N = c(nrow(vars$SideEffect),nrow(vars$GreaterGood))

    return(list(Response = factor(c(vars$SideEffect$haus1.1,vars$GreaterGood$haus2.1),levels=c(1,2),labels=vars$labels$Response),
                Condition = factor(c(rep(1,N[1]),rep(2,N[2])),levels=c(1,2),labels=vars$labels$Condition),
                N = N)
    )
}

varfun.Anderson.1 <- function(vars=ML2.sr){
    # Analysis plan. Following the original authors, the three dependent measures will be
    # standardized and averaged into a single index of subjective well-being. The mean difference in
    # subjective well-being between high and low-sociometric status conditions will be tested with an
    # independent-samples t-test. All participants with data will be included in the analysis.

    # and1.3=Satisfaction With Life Scale (SWLS, 5 items, Low SocioMetricStatus condition), higher numbers=higher satisfaction; and1.4=Positive And Negative Affect Scale (PANAS, Low SocioMetricStatus condition).
    # Positive items are 1,4,5,8,9,12,14,17,18,19. Negative items: 2,3,6,7,10,11,13,15,16,20;
    # Alert: recode responses to negative items before averaging.

    # and2.3=Satisfaction With Life Scale (SWLS, 5 items, High SocioMetricStatus condition), higher numbers=higher satisfaction; and2.4=Positive And Negative
    # Affect Scale (PANAS, High SocioMetricStatus condition). Positive items are 1,4,5,8,9,12,14,17,18,19.
    # Negative items: 2,3,6,7,10,11,13,15,16,20;
    # Alert: recode responses to negative items before averaging.

    #list(Low=c("and1.3_1", "and1.3_2", "and1.3_3", "and1.3_4", "and1.3_5", "and1.4_1", "and1.4_2", "and1.4_3", "and1.4_4", "and1.4_5", "and1.4_6", "and1.4_7", "and1.4_8", "and1.4_9", "and1.4_10", "and1.4_11", "and1.4_12", "and1.4_13", "and1.4_14", "and1.4_15", "and1.4_16", "and1.4_17", "and1.4_18", "and1.4_19", "and1.4_20"),
    # High=c("and2.3_1", "and2.3_2", "and2.3_3", "and2.3_4", "and2.3_5", "and2.4_1", "and2.4_2", "and2.4_3", "and2.4_4", "and2.4_5", "and2.4_6", "and2.4_7", "and2.4_8", "and2.4_9", "and2.4_10", "and2.4_11", "and2.4_12", "and2.4_13", "and2.4_14", "and2.4_15", "and2.4_16", "and2.4_17", "and2.4_18", "and2.4_19", "and2.4_20"))

    #
    # we computed a single subjective well-being score (SWB) by
    # standardising
    # then adding positive affect (alpha.0.84) and
    # life-satisfaction (alpha.0.84) and
    # subtracting the standardised negative affect variable (alpha.88).

    # Neagative
    PANASlowNA   <- select(vars$Low, one_of(unlist(strsplit(paste0("and1.4_",c(2,3,6,7,10,11,13,15,16,20),sep="|"),"|",fixed=T))))
    # Negative Recode
    PANASlowNAr  <- 6-PANASlowNA
    # Positive
    PANASlowPA   <- select(vars$Low, one_of(unlist(strsplit(paste0("and1.4_",c(1,4,5,8,9,12,14,17,18,19),sep="|"),"|",fixed=T))))
    # SWSL
    SWLSlow  <-  select(vars$Low, one_of(c("and1.3_1", "and1.3_2", "and1.3_3", "and1.3_4", "and1.3_5")))

    # Negative
    PANAShighNA  <- select(vars$High, one_of(unlist(strsplit(paste0("and2.4_",c(2,3,6,7,10,11,13,15,16,20),sep="|"),"|",fixed=T))))
    # Negative Recode
    PANAShighNAr <- 6-PANAShighNA
    # Positive
    PANAShighPA  <- select(vars$High, one_of(unlist(strsplit(paste0("and2.4_",c(1,4,5,8,9,12,14,17,18,19),sep="|"),"|",fixed=T))))
    # SWSL
    SWLShigh     <- select(vars$High,one_of(c("and2.3_1", "and2.3_2", "and2.3_3", "and2.3_4", "and2.3_5")))

    #   fiveSWLShigh <- scale.R(SWLShigh,lo=1,hi=5)
    #   fiveSWLSlow  <- scale.R(SWLSlow,lo=1,hi=5)

    # Descriptives
    MhighNA  <- rowMeans(scale(PANAShighNA))
    #MhighNAr <- rowMeans(scale(PANAShighNAr))
    MhighPA  <- rowMeans(scale(PANAShighPA))
    MhighSW  <- rowMeans(scale(SWLShigh))

    #SEShiD<- ldply(list(hiNA=MhighNA,hiNArec=MhighNAr,hiPA=MhighPA,highSW=MhighSW),function(m) data.frame(mean=mean(m),SD=sd(m),min=min(m),max=max(m)))

    MlowNA  <- rowMeans(scale(PANASlowNA))
    #MlowNAr <- rowMeans(scale(PANASlowNAr))
    MlowPA  <- rowMeans(scale(PANASlowPA))
    MlowSW  <- rowMeans(scale(SWLSlow))

    #  SESloD <- ldply(list(lowNA=MlowNA,lowNArec=MlowNAr,lowPA=MlowPA,lowSW=MlowSW),function(m) data.frame(mean=mean(m),SD=sd(m),min=min(m),max=max(m)))

    # GrandMean <- mean(c(rowMeans(cbind(PANASlowNAr,PANASlowPA,SWLSlow)),rowMeans(cbind(PANAShighNAr,PANAShighPA,SWLShigh))))

    return(list(Low  = (PANASlowPA+MlowSW)-MlowNA,
                High = (PANAShighPA+MhighSW)-MhighNA,
                #     Low  = rowMeans(cbind(sevenPANASlow,SWLSlow)),
                #               High = rowMeans(cbind(sevenPANAShigh,SWLShigh)),
                #     Low  = rowMeans(cbind(PANASlowNAr,PANASlowPA,SWLSlow))-GrandMean,
                #               High = rowMeans(cbind(PANAShighNAr,PANAShighPA,SWLShigh))-GrandMean,
                N    = vars$N))
    #               GrandMean = GrandMean,
    #               Data = list(LowSES =data.frame(PANAS=cbind(PANASlowNAr,PANASlowPA),SWLS=SWLSlow),
    #                           HighSES=data.frame(PANAS=cbind(PANAShighNAr,PANAShighPA),SWLS=SWLShigh)
    #               ),
    #               Descriptives = list(LowSESdesc  <- SESloD,
    #                                   HighSESdesc <- SEShiD)

}

varfun.Ross.1 <- function(vars=ML2.sr, DEMO=FALSE){
    #   Analysis plan. An independent samples t-test will be conducted with participants’ choice
    # (sign release/refuse to sign) as the IV and participant estimate of % of peers who would sign the
    # release as the DV. Note that participants self-select whether to sign or refuse to sign the release,
    # so it is not random assignment to levels of the independent variable. Participants will be included
    # in the analysis if they respond to all three questions and their estimate for the DV (e.g., “what
    # percent of your peers would sign the release”) falls between 0-100
    #ross.s1.1 = percentage of peers;
    #ross.s1.2 = you; values: 1=sign; 2=refuse
    # list(Peers='ross.s1.1',You='ross.s1.2')

    if(DEMO==TRUE){

        warning('Ross.1: DEMO=TRUE',call. = FALSE)
        return(list(Peers = rnorm(50),
                    You   = c(rep(1,25),rep(2,25)),
                    N     = c(50,50)))

    } else {

        return(list(Peers  = vars$Peers[[1]],
                    You    = factor(vars$Peers[[2]],levels=c(1,2),labels=vars$labels$You),
                    N      = c(sum(vars$Peers[[2]]==1),sum(vars$Peers[[2]]==2)))
        )}
}

varfun.Ross.2 <- function(vars=ML2.sr){
    # Analysis plan. An independent samples t-test will be conducted with participants’ choice
    # (pay the fine/appear in court) as the IV and participant estimate of percent of peers who would
    # pay the fine as the DV. Note that participants self-select whether to pay the fine or appear in
    # court, so it is not random assignment to levels of the independent variable. Participants will be
    # included in the analysis if they respond to all three questions and their estimate for the DV (e.g.,
    # “what percent of your peers would pay the fine by mail”) falls between 0-100

    #ross.s2.1=percentage of peers;
    #ross.s2.2=you; values: 1=Pay; 2=Appear in court
    #  list(Peers='ross.s2.1_1_TEXT',You='ross.s2.2')

    return(list(Peers  = vars$Peers[[1]],
                You    = factor(vars$Peers[[2]],levels=c(1,2),labels=vars$labels$You),
                N      = c(sum(vars$Peers[[2]]==1), sum(vars$Peers[[2]]==2)))
    )
}

varfun.Giessner.1 <- function(vars=ML2.sr){
    # Analysis plan. Responses to the five dependent measures will be averaged and an
    # independent samples t-test will compare mean power rating between the 2 cm and 7 cm
    # conditions. All participants who complete at least one item of the dependent measure will be
    # included in the analysis.

    # geis.1.1=long line condition;
    # geis.2.1=short line condition;
    # geis.dv_1=dominant;
    # geis.dv_2=strong;
    # geis.dv_3=self-confident;
    # geis.dv_4=control;
    # geis.dv_5=status;
    # For all dvs, higher numbers=higher power.

    # list(Line=c("geis.1.1","geis.2.1"),Power=c("geis.dv_1","geis.dv_2","geis.dv_3","geis.dv_4","geis.dv_5"))

    #   id.ValidL <- ML2.id$&ML2.id$You&(ML2.df[ML2.in$study.vars$You]==1)
    #   id.ValidS <- ML2.id$Peers&ML2.id$You&(ML2.df[ML2.in$study.vars$You]==2)
    Long <- vars$Long  %>% filter(!is.na(vars$Long[1])  & rowSums(!is.na(vars$Long[-1, ]))>1)
    Short<- vars$Short %>% filter(!is.na(vars$Short[1]) & rowSums(!is.na(vars$Short[-1, ]))>1)

    return(list(Long   = rowMeans(Long[-1]),
                Short  = rowMeans(Short[-1]),
                N      = c(nrow(Long),nrow(Short)))
    )
}

varfun.Tversky.1 <- function(vars=ML2.sr){
    # Tversky, A., Kahneman, D. (1981). The framing of decisions and the psychology of choice. Science, 211, 453-458.
    # tver1.1=choice ($250 wall hanging condition, yes=1, no=2);
    # tver2.1=choice ($30 wall hanging cond, yes=1, no=2).

    #   Materials and procedure. Participants will receive one of two scenarios from the original with dollar amounts approximately adjusted for inflation and the consumer items being replaced with a ceramic vase and a wall hanging.

    # tver1.1 Imagine that you are about to purchase a ceramic vase for $30, and a wall hanging for $250. The salesman informs you that the wall hanging you wish to buy is on sale for $240 at the other branch of the store, located 20 minutes drive away. Would you make the trip to the other store?
    # m  Yes, I would go to the other branch. (1)
    # m  No, I would not go to the other branch. (2)

    # tver2.1 Imagine that you are about to purchase a ceramic vase for $250, and a wall hanging for $30. The salesman informs you that the wall hanging you wish to buy is on sale for $20 at the other branch of the store, located 20 minutes drive away. Would you make the trip to the other store?
    # m  Yes, I would go to the other branch. (1)
    # m  No, I would not go to the other branch. (2)

    # Materials here: https://osf.io/8t9ha/. Test the study here: https://ufl.qualtrics.com/SE/?SID=SV_aW8rIyGPNQ2Wwh7.

    # Analysis plan. A two-way contingency table will be built with Price condition ($20 vs. $240) and Response (Yes vs. No) as factors. The critical replication hypothesis will be given by a χ2 test and the effect size by an odds ratio. All participants with valid responses will be included in analysis.

    return(list(Condition = factor(c(vars$Vase[[1]],vars$Wall[[1]]),levels=c(1,2),labels=vars$labels$Condition),
                Response  = factor(c(rep(1,vars$N[1]),rep(2,vars$N[2])),levels=c(1,2),labels=vars$labels$Response),
                N      = vars$N)
    )
}


varfun.Hauser.2 <- function(vars=ML2.sr){
    # Hauser, M., Cushman, F., Young, L., Kang-Xing Jin, R., & Mikhail, J. (2007). A dissociation between moral judgments and justifications. Mind & Language, 22, 1-21.

    #   hauser3.1=morality judgment (greater good condition);
    #   hauser4.1=morality judgment (foreseen side-effect condition; for both, yes=1, no=2.)
    #   haus3.2 and haus4.2=previous experience (yes=1, no=2);
    #   haus3.1t_3=timing (greater good);
    #   haus4.1t_3=timing (side effect).


    # Participants will be excluded from all analyses if they take fewer than four seconds to read and respond to either of the target scenarios. For the key confirmatory test comparing with the original effect size, we will include only participants that indicate having no prior experience with the task. Prior exposure will be investigated as a moderator for the other analyses

    SideEffect  <- filter(vars$SideEffect, vars$SideEffect[vars$labels$Experience[1]] !=1 & vars$SideEffect[vars$labels$Timing[1]]>4)
    GreaterGood <- filter(vars$GreaterGood,vars$GreaterGood[vars$labels$Experience[2]]!=1 & vars$GreaterGood[vars$labels$Timing[2]]>4)
    N <- c(nrow(SideEffect),nrow(GreaterGood))

    return(list(Response  = factor(c(SideEffect$haus3.1,GreaterGood$haus4.1),levels=c(1,2),labels=vars$labels$Response),
                Condition = factor(c(rep(1,N[1]),rep(2,N[2])),levels=c(1,2),labels=vars$labels$Condition),
                N = N)
    )
}

varfun.Risen.1 <- function(vars=ML2.sr){
    # Risen, J. L., & Gilovich, T. (2008). Why people are reluctant to tempt fate. Journal of Personality and Social Psychology, 95, 293.

    # rise1.3=likelihood that the professor will call on you (unprepared condition);
    # rise2.3=likelihood that the professor will call on you (prepared condition);

    # for both, higher numbers=higher likelihood	 All participants that answer the dependent measure will be included in analysis.
    # The primary confirmatory test for comparing the original and replication effect size will be based on only the samples using undergraduate students.

    #   Analysis plan. The two groups will be compared with an independent samples t-test. All participants that answer the dependent measure will be included in analysis.  The primary confirmatory test for comparing the original and replication effect size will be based on only the samples using undergraduate students. We will examine gender as a possible moderator of the effect in a supplemental, exploratory analysis.

    #Variable = "ex.subjp" which asked if participants were recruited through a university subject pool. 1 = yes, 2 = no.

    return(list(Unprepared  = as.numeric(vars$Unprepared[[1]]),
                Prepared    = as.numeric(vars$Prepared[[1]]),
                N           = vars$N))
}


varfun.Risen.2 <- function(vars=ML2.sr){
    # Risen, J. L., & Gilovich, T. (2008). Why people are reluctant to tempt fate. Journal of Personality and Social Psychology, 95, 293.

    # rise1.3=likelihood that the professor will call on you (unprepared condition);
    # rise2.3=likelihood that the professor will call on you (prepared condition);

    # for both, higher numbers=higher likelihood   All participants that answer the dependent measure will be included in analysis.
    # The primary confirmatory test for comparing the original and replication effect size will be based on only the samples using undergraduate students.

    #   Analysis plan. The two groups will be compared with an independent samples t-test. All participants that answer the dependent measure will be included in analysis.  The primary confirmatory test for comparing the original and replication effect size will be based on only the samples using undergraduate students. We will examine gender as a possible moderator of the effect in a supplemental, exploratory analysis.

    #use variable "sex" 1 = male, 2 = female. (3 = other and 4 = prefer not to answer, which I think we probably do not use for this analysis)

    return(list(Likelihood = c(vars$Unprepared$rise1.3,vars$Prepared$rise2.3),
                Condition  = factor(c(rep(1,vars$N[1]),rep(2,vars$N[2])),levels=c(1,2),labels=vars$labels$Condition),
                Gender     = factor(c(vars$Unprepared$sex,vars$Prepared$sex)),
                N          = vars$N)
    )
}

varfun.Savani.1 <- function(vars=ML2.sr){
    # Savani, K., Markus, H. R., Naidu, N. V. R., Kumar, S., & Berlia, N. (2010). What counts as a choice? US Americans are more likely than Indians to construe actions as choices. Psychological Science, 21, 391-398.

    # sava1.N=interpersonal actions;
    # sava2.N=personal actions;

    # 'sava1.4', 'sava1.5', 'sava1.9', 'sava1.10', 'sava1.15', 'sava1.16' , 'sava1.21', 'sava1.22', 'sava1.27', 'sava1.28', 'sava1.33', 'sava1.34',  'sava1.38', 'sava1.39', 'sava1.43', 'sava1.44'

    # sava1.4=choice (buy a gift; 1=choice; 2=no choice);
    # sava1.5=importance (buy a gift);
    # sava1.9=choice (take a friend at the restaurant; 1=choice; 2=no choice);
    # sava1.10=importance (restaurant);
    # sava1.15=choice (trip; 1=choice, 2 and 3 = no choice);
    # sava1.16=importance (trip);
    # sava1.21=choice (dinner; 1=choice, 2 and 3 = no choice);
    # sava1.22=importance (dinner);
    # sava1.27=choice (errand; 1=choice, 2 and 3 = no choice);
    # sava1.28=importance (errand);
    # sava1.33=choice (help, 1=choice, 2 & 3 = no choice);
    # sava1.34=importance (help);
    # sava1.38=choice (advice, 1=choice, 2 & 3 = no choice);
    # sava1.39=importance (advice);
    # sava1.43=choice (friends, 1=choice, 2 & 3 = no choice);
    # sava1.44=importance (friends);

    # sava2.4=choice (buy for yourself; 1=choice; 2=no choice);
    # sava2.5=importance (buy for yourself);
    # sava2.9=choice (at the restaurant by yourself; 1=choice; 2=no choice);
    # sava2.10=importance (restaurant by yourself);
    # sava2.15=choice (trip alone; 1=choice, 2 and 3 = no choice);
    # sava2.16=importance (trip alone);
    # sava2.21=choice (out for dinner; 1=choice, 2 and 3 = no choice);
    # sava2.22=importance (out for dinner);
    # sava2.27=choice (errand for yourself; 1=choice, 2 and 3 = no choice);
    # sava2.28=importance (errand for yourself);
    # sava2.33=choice (ask for help, 1=choice, 2 & 3 = no choice);
    # sava2.34=importance (ask for help);
    # sava2.38=choice (take a course, 1=choice, 2 & 3 = no choice);
    # sava2.39=importance (take a course);
    # sava2.43=choice (friends, 1=choice, 2 & 3 = no choice);
    # sava2.44=importance (friends);

    # For all importance items: higher numbers=higher importance

    # we will only include university data collections in the primary confirmatory analysis to be compared with the original effect sizes.

    # Data for all participants will be included to examine variability across sample and setting.
    # However, participants must respond to all choice and importance of choice questions to be included in the analysis.

    #   Analysis plan. We will conduct a hierarchical logistic regression analysis with Choice (binary) as the dependent variable, the Importance of decision (ordered categorical) as a trial-level covariate nested within participants, and Condition (categorical) as a participant-level factor indicating whether a participant was in the personal or interpersonal condition.  The effect of interest will be the odds of construing an action as a choice, depending on the condition a participant was in, controlling for the reported importance of the action.

    # Because some survey questions may be less suitable for non-student samples, we will only include university data collections in the primary confirmatory analysis to be compared with the original effect sizes. Data for all participants will be included to examine variability across sample and setting. However, participants must respond to all choice and importance of choice questions to be included in the analysis. The target effect size for replication will be the the results obtained for participants from labs in India, to compare to the effect found in the Indian sample in the original (Indian participants were more likely to construe interpersonal actions as choices than personal actions). Although we only have few labs from India, we are making extra efforts to recruit many participants in those labs. We anticipate that this effect will vary by sample, particularly in line with the original demonstration of cultural differences.
    #
    in.IT(c('lme4','lmerTest'))

    choice.Int     <- c('sava1.4', 'sava1.9', 'sava1.15', 'sava1.21', 'sava1.27', 'sava1.33', 'sava1.38', 'sava1.43')
    importance.Int <- c('sava1.5', 'sava1.10', 'sava1.16', 'sava1.22', 'sava1.28', 'sava1.34', 'sava1.39', 'sava1.44')

    choice.Pers     <- c('sava2.4', 'sava2.9', 'sava2.15', 'sava2.21', 'sava2.27', 'sava2.33', 'sava2.38', 'sava2.43')
    importance.Pers <- c('sava2.5','sava2.10', 'sava2.16', 'sava2.22', 'sava2.28','sava2.34', 'sava2.39', 'sava2.44')

    Condition    <- factor(c(rep(1,vars$N[1]),rep(2,vars$N[2])),levels=c(1,2),labels=vars$labels$Condition)
    id           <- seq_along(Condition)

    Response             <- rbind(dplyr::select(vars$Interpersonal,Choice=one_of(choice.Int)), dplyr::select(vars$Personal,Choice=one_of(choice.Pers)))
    Response[Response>1] <- 2
    Response$uID         <- id
    Response$Condition   <- Condition
    Response             <- melt(Response,id=c('uID','Condition'),variable.name='trialID',value.name='Response')
    Response$Response    <- factor(Response$Response,levels=c(1,2),labels=vars$labels$Response)

    Importance           <- rbind(dplyr::select(vars$Interpersonal,Importance=one_of(importance.Int)), dplyr::select(vars$Personal,Importance=one_of(importance.Pers)))
    Importance$uID       <- id
    Importance$Condition <- Condition
    Importance           <- melt(Importance,id=c('uID','Condition'),variable.name='VarLabel',value.name='Importance')
    # Probably superfluous
    matchID <- Response$uID==Importance$uID
    Response$Importance[matchID]           <- scale(Importance$Importance[matchID],scale=F)

    #     m0<-   glmer(Response ~ Importance*Condition + (1|uID/trialID),family = binomial('logit'),data=df)
    #     exp(fixef(m0)[4])

    return(list(df = tbl_df(Response),
                N  = vars$N))
}

varfun.Savani.2 <- function(vars=ML2.sr){
    # Savani, K., Markus, H. R., Naidu, N. V. R., Kumar, S., & Berlia, N. (2010). What counts as a choice? US Americans are more likely than Indians to construe actions as choices. Psychological Science, 21, 391-398.

    # Exact copy of Savani.1 for secondary analysis:
    # we will only include university data collections in the primary confirmatory analysis to be compared with the original effect sizes.

    # Because some survey questions may be less suitable for non-student samples, we will only include university data collections in the primary confirmatory analysis to be compared with the original effect sizes. Data for all participants will be included to examine variability across sample and setting. However, participants must respond to all choice and importance of choice questions to be included in the analysis. The target effect size for replication will be the the results obtained for participants from labs in India, to compare to the effect found in the Indian sample in the original (Indian participants were more likely to construe interpersonal actions as choices than personal actions). Although we only have few labs from India, we are making extra efforts to recruit many participants in those labs. We anticipate that this effect will vary by sample, particularly in line with the original demonstration of cultural differences.

    in.IT(c('lme4','lmerTest'))

    choice.Int     <- c('sava1.4', 'sava1.9', 'sava1.15', 'sava1.21', 'sava1.27', 'sava1.33', 'sava1.38', 'sava1.43')
    importance.Int <- c('sava1.5', 'sava1.10', 'sava1.16', 'sava1.22', 'sava1.28', 'sava1.34', 'sava1.39', 'sava1.44')

    choice.Pers     <- c('sava2.4', 'sava2.9', 'sava2.15', 'sava2.21', 'sava2.27', 'sava2.33', 'sava2.38', 'sava2.43')
    importance.Pers <- c('sava2.5','sava2.10', 'sava2.16', 'sava2.22', 'sava2.28','sava2.34', 'sava2.39', 'sava2.44')

    Condition    <- factor(c(rep(1,vars$N[1]),rep(2,vars$N[2])),levels=c(1,2),labels=vars$labels$Condition)
    id           <- seq_along(Condition)

    Response             <- rbind(dplyr::select(vars$Interpersonal,Choice=one_of(choice.Int)), dplyr::select(vars$Personal,Choice=one_of(choice.Pers)))
    Response[Response>1] <- 2
    Response$uID         <- id
    Response$Condition   <- Condition
    Response             <- melt(Response,id=c('uID','Condition'),variable.name='trialID',value.name='Response')
    Response$Response    <- factor(Response$Response,levels=c(1,2),labels=vars$labels$Response)

    Importance           <- rbind(dplyr::select(vars$Interpersonal,Importance=one_of(importance.Int)), dplyr::select(vars$Personal,Importance=one_of(importance.Pers)))
    Importance$uID       <- id
    Importance$Condition <- Condition
    Importance           <- melt(Importance,id=c('uID','Condition'),variable.name='VarLabel',value.name='Importance')
    # Probably superfluous
    matchID <- Response$uID==Importance$uID
    Response$Importance[matchID]           <- scale(Importance$Importance[matchID],scale=F)

    return(list(df = tbl_df(Response),
                N  = vars$N))
}


varfun.Norenzayan.1 <- function(vars=ML2.sr){
    # Norenzayan, A., Smith, E. E., Kim, B. J., & Nisbett, R. E. (2002). Cultural preferences for formal versus intuitive reasoning. Cognitive Science, 26, 653-684.
    # "nore1.1 TO nore1.20 provide choices (""belong to"" condition);
    # nore2.1 to nore2.20 provide choices (""similar to"" condition).
    # [WAIT FOR APPROVAL OF SCORING KEYS]"

    # # Analysis plan. We will compute for each subject the percentage of rule-based responses and test whether the mean of the two experimental groups (belong vs similar) on this DV is equal with a t-test for independent samples. The effect size will be given by a standardized mean difference. All participants with data will be included in analysis.
    # # For additional analysis, a few items about cultural origins of participants and their parents are present in the individual differences assessment.  These could be particularly useful for follow-up moderator analysis.

    return(list(Belong  = mean(as.numeric(vars$Belong[[1]],na.rm=T)),
                Similar = mean(as.numeric(vars$Similar[[1]],na.rm=T)),
                N       = vars$N))
}

varfun.Hsee.1 <- function(vars=ML2.sr){
    # Hsee, C. K. (1998). Less is better: When low-value options are valued more highly than high-value options. Journal of Behavioral Decision Making, 11, 107-121.

    #hsee1.1=generosity ($90 scarf condition);
    #hsee2.1=generosity ($110 coat condition);
    #for both, higher numbers=higher generosity

    #   Analysis plan. The two conditions will be compared with an independent samples t-test
    # with rated generosity of gift giver as the dependent variable. All participants with data will be
    # included in the analysis
    return(list(Scarf  = as.numeric(vars$Scarf[[1]]),
                Coat   = as.numeric(vars$Coat[[1]]),
                N      = vars$N))
}

varfun.Gray.1 <- function(vars=ML2.sr){
    # Gray, K., & Wegner, D. M. (2009). Moral typecasting: divergent perceptions of moral agents and moral patients. Journal of Personality and Social Psychology, 96, 505.
    # Adult harms baby scenario; gray1.2=responsibility (adult);  gray1.4=pain (baby).
    # Baby harms adult scenario; gray2.2=responsibility (baby);  gray2.4=pain (adult)

    #   Analysis plan. We will compare the means on perceived responsibility between
    # conditions with an independent samples t-test for the responsibility item. The intentionality item
    # will be analyzed the same way as a secondary analysis. All participants with data will be
    # included in analysis.

    return(list(adultHbaby = as.numeric(vars$adultHbaby[[1]]),
                babyHadult = as.numeric(vars$babyHadult[[1]]),
                N     = vars$N))
}


varfun.Gray.2 <- function(vars=ML2.sr){
    # Gray, K., & Wegner, D. M. (2009). Moral typecasting: divergent perceptions of moral agents and moral patients. Journal of Personality and Social Psychology, 96, 505.
    # Baby harms adult scenario; gray1.3=intentionality (adult); gray1.4=pain (baby).
    # Adult harms baby scenario; gray2.3=intentionality (baby); gray2.4=pain (adult)

    #   Analysis plan. We will compare the means on perceived responsibility between
    # conditions with an independent samples t-test for the responsibility item.

    # The intentionality item will be analyzed the same way as a secondary analysis. All participants with data will be
    # included in analysis.

    return(list(adultHbaby = as.numeric(vars$adultHbaby[[1]]),
                babyHadult = as.numeric(vars$babyHadult[[1]]),
                N     = vars$N))
}


varfun.Zhong.1 <- function(vars=ML2.sr){
    # Zhong, C. B., & Liljenquist, K. (2006). Washing away your sins: Threatened morality and physical cleansing. Science, 313, 1451–1452.

    # zhon1.1= unethical condition;
    # zhon2.1=ethical condition;
    # zhon.dv.1_1 TO zhon.dv.1_10=desirability of products (both conditions); higher numbers=higher desirability; # Products list:
    # Post-it notes (zhon.dv.1_1),
    # Dove shower soap (zhon.dv.1_2),
    # Crest toothpaste (zhon.dv.1_3),
    # Nantucket Nectars juice (zhon.dv.1_4),
    # Energizer batteries (zhon.dv.1_5),
    # Sony cd cases (zhon.dv.1_6),
    # Windex glass cleaner (zhon.dv.1_7),
    # Lysol countertop disinfectant (zhon.dv.1_8),
    # Snickers candy bar (zhon.dv.1_9),
    # Tide laundry detergent (zhon.dv.1_10)

    # Participants who copy less than half the target article will be excluded from analysis

    #   Analysis plan. The key factor of interest is whether condition affects ratings of the
    # cleaning products, so ratings of the five cleaning products will be averaged and compared
    # between the two conditions (ethical or unethical) with an independent samples t-test. A second
    # comparison is whether there is a condition difference between ratings of the other products. The
    # theoretical expectation is that this difference will be weak or near zero. This will be examined as
    # a secondary analysis as a 2 (ethical-unethical) x 2 (cleaning-other products) mixed model
    # ANOVA, and a follow-up independent samples t-test comparing ratings of other products
    # between the ethical and unethical conditions.3

    # Participants who copy less than half the target article will be excluded from analysis.
    return(list(Ethical   = rowMeans(vars$Ethical[,2:11]),
                Unethical = rowMeans(vars$Unethical[,2:11]),
                N         = vars$N))
}


varfun.Zhong.2 <- function(vars=ML2.sr){
    # Zhong, C. B., & Liljenquist, K. (2006). Washing away your sins: Threatened morality and physical cleansing. Science, 313, 1451–1452.

    # zhon1.1= unethical condition;
    # zhon2.1=ethical condition;

    # zhon.dv.1_1 TO zhon.dv.1_10=desirability of products (both conditions); higher numbers=higher desirability; # Products list:

    # Clean:
    # Dove shower soap (zhon.dv.1_2),
    # Crest toothpaste (zhon.dv.1_3),
    # Windex glass cleaner (zhon.dv.1_7),
    # Lysol countertop disinfectant (zhon.dv.1_8),
    # Tide laundry detergent (zhon.dv.1_10)

    # Not-clean:
    # Post-it notes (zhon.dv.1_1),
    # Nantucket Nectars juice (zhon.dv.1_4),
    # Energizer batteries (zhon.dv.1_5),
    # Sony cd cases (zhon.dv.1_6),
    # Snickers candy bar (zhon.dv.1_9),

    # Participants who copy less than half the target article will be excluded from analysis

    # Analysis plan. The key factor of interest is whether condition affects ratings of the
    # cleaning products, so ratings of the five cleaning products will be averaged and compared
    # between the two conditions (ethical or unethical) with an independent samples t-test. A second
    # comparison is whether there is a condition difference between ratings of the other products. The
    # theoretical expectation is that this difference will be weak or near zero. This will be examined as


    # a secondary analysis as a 2 (ethical-unethical) x 2 (cleaning-other products) mixed model
    # ANOVA, and a follow-up independent samples t-test comparing ratings of other products
    # between the ethical and unethical conditions.3
    # Participants who copy less than half the target article will be excluded from analysis.
    in.IT(c('lme4','lmerTest'))

    idClean <- c(2,3,7,8,10)
    idOther <- c(1,4,5,6,9)

    return(list(Response  = c(rowMeans(vars$Ethical[,idClean+1]), rowMeans(vars$Ethical[,idOther+1]), rowMeans(vars$Unethical[,idClean+1]), rowMeans(vars$Unethical[,idOther+1])),
                Condition = factor(c(vars$Ethical[,1],vars$Ethical[,1], vars$Unethical[,1],vars$Unethical[,1]),levels=c(2,1),labels=vars$labels$Condition),
                Product   = factor(c(rep(1,times=vars$N[1]),rep(2,times=vars$N[1]),rep(1,times=vars$N[2]),rep(2,times=vars$N[2])),levels=c(1,2),labels=vars$labels$Product),
                uID       = c(seq_along(vars$Ethical[,1]),seq_along(vars$Ethical[,1]),max(seq_along(vars$Ethical[,1]))+seq_along(vars$Unethical[,1]),max(seq_along(vars$Ethical[,1]))+seq_along(vars$Unethical[,1])),
                N         = vars$N))
}


varfun.Schwarz.1 <- function(vars=ML2.sr){
    # Schwarz, N., Strack, F., & Mai, H. P. (1991). Assimilation and contrast effects in part-whole question sequences: A conversational logic analysis. Public Opinion Quarterly, 55, 3-23.

    # schw1.1=life sat (first);
    # schw1.2=partner satisfaction (second);
    # schw2.1=partner satisfaction (first);
    # schw2.2=life sat (second).
    # for all, higher numbers=higher satisfaction

    #   Analysis plan. We will compute the correlation between responses to the general and
    # specific question in each item order condition, and then compare the correlations using the Fisher
    # r-to-z transformation. Participants with valid responses to both items will be included in the
    # analysis.
    return(list(Fisher = list(r1=cor(vars$SpecificFirst[,1],vars$SpecificFirst[,2],use="pairwise.complete.obs"),
                              r2=cor(vars$GlobalFirst[,1],vars$GlobalFirst[,2],use="pairwise.complete.obs"),
                              n1=vars$N[1],
                              n2=vars$N[2]),
                N = vars$N))

}


varfun.Shafir.1 <- function(vars=ML2.sr){
    # Shafir, E. (1993). Choosing versus rejecting: Why some options are both better and worse than others. Memory & Cognition, 21, 546-556.

    # shaf1.1=choice (award condition; Parent A=1, Parent B=2);
    # shaf2.1=choice (deny condition; Parent A=1, Parent B=2)
    #
    # Analysis plan. The proportion of participants awarding or denying custody for parent B
    # will be summed from both groups and tested against 100% with a Z test. The effect size will be
    # computed estimating a logistic regression model on the 2X2 table and then taking the
    # exponentiation of the unstandardized beta parameter of the main effect of parent B, which can be
    # interpreted as an odds ratio. All participants with data will be included in analysis.



}

varfun.Zaval.1 <- function(vars=ML2.sr){
    # Zaval, L., Keenan, E. A., Johnson, E. J., & Weber, E. U. (2014). How warm days increase belief in global warming. Nature Climate Change, 4, 143-147.

    #zav1.1 TO zav1.13 provide COLD primes.
    #zav2.1 TO zav2.13 provide HEAT primes.
    #zav.dv.2=belief;
    #zav.dv.3=concern;
    #higher numbers=higher belief/concern.

    #the direct replication test will use only English language sites, and - like all other effects - all samples and settings with data will be included in analyses examining heterogeneity to see if factors, like translation, have an impact on effect estimates

    #   Analysis plan. Mean differences in belief and concern about global warming between
    # heat and cold-priming conditions will be evaluated with an independent samples t-test. The
    # scrambled sentence task could introduce unanticipated variation across translations. As such, the
    # direct replication test will use only English language sites, and - like all other effects - all
    # samples and settings with data will be included in analyses examining heterogeneity to see if
    # factors, like translation, have an impact on effect estimates.

    return(list(Cold = c(vars$Cold[,2]-vars$Cold[,3]),
                Heat = c(vars$Heat[,2]-vars$Heat[,3]),
                N    = vars$N))

}


varfun.Zaval.2 <- function(vars=ML2.sr){# ONLY ENGLISH
    # Zaval, L., Keenan, E. A., Johnson, E. J., & Weber, E. U. (2014). How warm days increase belief in global warming. Nature Climate Change, 4, 143-147.

    #zav1.1 TO zav1.13 provide COLD primes.
    #zav2.1 TO zav2.13 provide HEAT primes.
    #zav.dv.2=belief;
    #zav.dv.3=concern;
    #higher numbers=higher belief/concern.

    #the direct replication test will use only English language sites, and - like all other effects - all samples and settings with data will be included in analyses examining heterogeneity to see if factors, like translation, have an impact on effect estimates

    #   Analysis plan. Mean differences in belief and concern about global warming between
    # heat and cold-priming conditions will be evaluated with an independent samples t-test. The
    # scrambled sentence task could introduce unanticipated variation across translations. As such, the
    # direct replication test will use only English language sites, and - like all other effects - all
    # samples and settings with data will be included in analyses examining heterogeneity to see if
    # factors, like translation, have an impact on effect estimates.

    return(list(Cold = c(vars$Cold[,2]-vars$Cold[,3]),
                Heat = c(vars$Heat[,2]-vars$Heat[,3]),
                N    = vars$N))

}

varfun.Knobe.1 <- function(vars=ML2.sr){
    # Knobe, J. (2003). Intentional action and side effects in ordinary language. Analysis, 63, 190-193.

    #knob1.3=intentionality (help condition);
    #knob2.3=intentionality (harm condition);
    #for both, higher numbers=higher intentionality

    #   Analysis plan. Ratings of intentionality in the harm and help conditions will be
    # compared using an independent-samples t-test. This is the focal test for the direct replication.
    # Blame and praise ratings will also be collected but are secondary analyses. All participants with
    # data will be included in analysis.

    return(list(Help = vars$Help,
                Harm = vars$Harm,
                N    = vars$N))
}

varfun.Gati.1 <- function(vars=ML2.sr){
    # Tversky, A., & Gati, I. (1978). Studies of similarity. Cognition and categorization, 1, 79-98.
    # see Gati sheet

    #   Analysis plan. We will perform three analyses on the data.

    # For the primary analysis, we
    # will analyze the data through a general linear mixed model with a random effect for the item pair
    # nested within subject, and a fixed factor ‘order’ representing the order of the pair (prominent first
    # vs. prominent second). Fitting this model will allow evaluation of both effects. If the intercept is
    # significantly greater than 0, this would confirm the finding that at the participant level, if there is
    # an effect for the factor ‘order’ the pairs where the prominent country appeared second will be
    # rated as more similar than when the prominent country appeared first. We will convert the Beta
    # provided by this intercept term into a Cohen’s d effect size.

    # Second, we will recreate the original analysis used to get a participant-level effect of
    # making similarity judgments where either the more or less prominent country comes first. We
    # will compute an asymmetry score for each subject, calculated as the average similarity for
    # comparisons where the prominent country appears second minus the average for the comparisons
    # where the prominent country appears first. Using a one-sample t-test, we will test this difference
    # score against zero (original d=.48).

    # Third, using a matched-pairs t-test, we will compare the
    # average score for each pair when it was prominent-first compared to prominent-second (original d = 1.31).

    # Because these latter two analyses do not account for the fact that the variance in ratings is
    # crossed between participants and pairs, they will be secondary and only used as a comparison for
    # the original analysis. All participants with data will be included in the analysis.
    # These analyses will be repeated for the differences conditions and reported as a separate
    # study. Because of the random assignment to similarity or difference conditions, each site will
    # have half as much data for its critical test as the other effects. This will likely increase the
    # standard error of its estimates by comparison.

    # PromFirst  gati2	4	6	7	8	9	10	13	15	20	21
    # PromSecond	gati2	2	3	5	11	12	14	16	17	18	19	22
    # PromFirst	gati1	2	3	5	11	12	14	16	17	18	19	22
    # PromSecond	gati1	4	6	7	8	9	10	13	15	20	21

    # CB1 <- list(P1st = c(2,	3,	5,	11,	12,	14,	16,	17,	18,	19,	22),
    #             P2nd = c(4,  6,	7,	8,	9,	10,	13,	15,	20,	21))
    #
    # CB1Diff <- llply(CB1,function(i) paste('gati1',i,sep="d."))
    # CB1Sim  <- llply(CB1,function(i) paste('gati1',i,sep="s."))
    #
    # CB1DiffV <- llply(CB1Diff,function(i) paste0("c('",paste(i,collapse="','"),"')"))
    # CB1SimV  <- llply(CB1Sim,function(i) paste0("c('",paste(i,collapse="','"),"')"))
    #
    # CB2 <- list(P1st = c(4,  6,  7,	8,	9,	10,	13,	15,	20,	21),
    #             P2nd = c(2,  3,  5,	11,	12,	14,	16,	17,	18,	19,	22))
    #
    # CB2Diff  <- llply(CB2,function(i) paste("gati2",i,sep="d."))
    # CB2Sim   <- llply(CB2,function(i) paste("gati2",i,sep="s."))
    #
    # CB2DiffV <- llply(CB2Diff,function(i) paste0("c('",paste(i,collapse="','"),"')"))
    # CB2SimV  <- llply(CB2Sim,function(i) paste0("c('",paste(i,collapse="','"),"')"))

    vars$PrimeFirstA$uID       <- seq_along(vars$PrimeFirstA[,1])
    vars$PrimeFirstA$Condition <- 1
    vars$PrimeFirstA$CounterB  <- 1
    vars$PrimeFirstB$uID       <- seq_along(vars$PrimeFirstB[,1])+max(vars$PrimeFirstA$uID)
    vars$PrimeFirstB$Condition <- 1
    vars$PrimeFirstB$CounterB  <- 2

    Similarity1 <- melt(vars$PrimeFirstA,id=c('uID','Condition','CounterB'),variable.name='itemID',value.name='Similarity')

    vars$PrimeSecondA$uID      <- seq_along(vars$PrimeSecondA[,1])+max(vars$PrimeFirstB$uID)
    vars$PrimeSecondA$Condition<- 2
    vars$PrimeSecondA$CounterB <- 1
    vars$PrimeSecondB$uID      <- seq_along(vars$PrimeSecondB[,1])+max(vars$PrimeSecondA$uID)
    vars$PrimeSecondB$Condition<- 2
    vars$PrimeSecondB$CounterB <- 2

    Similarity2 <- melt(vars$PrimeSecondB,id=c('uID','Condition','CounterB'),variable.name='itemID',value.name='Similarity')

    df <- as.data.frame(rbind(Similarity1,Similarity2))

    df$Condition <- factor(df$Condition,levels=c(1,2),labels=c('PrimeFirst','PrimeSecond'))
    df$CounterB  <- factor(df$CounterB,levels=c(1,2),labels=c('CntrBalanceA','CntrBalanceB'))
    itemID       <- unlist(strsplit(as.character(df$itemID),".", fixed=T)) #factor(df$itemID)
    df$itemID    <- as.numeric(itemID[grepl('(\\b(\\d)+\\b)',itemID)])

    vars$N <- c(vars$N[1]+vars$N[2],vars$N[3]+vars$N[4])

    # M0 <- lmer(Similarity ~ Condition + (1|uID) + (1|itemID),data=df)

    return(list(df = tbl_df(df),
                N = vars$N))

}
#
# Similar   <- paste0("c('",paste(colnames(ML2.S2)[grepl('(gati)(\\d)(s[.])(\\d)',colnames(ML2.S2))],collapse="','"),"')")
# Different <- paste0("c('",paste(colnames(ML2.S2)[grepl('(gati)(\\d)(d[.])(\\d)',colnames(ML2.S2))],collapse="','"),"')")


#
# (select(ML2.S2,starts_with('gati1s.')))[1:20,]
# (select(ML2.S2,starts_with('gati1d.')))[1:20,]
#
# (select(ML2.S2,starts_with('gati2s.')))[1:10,]
# (select(ML2.S2,starts_with('gati2d.')))[1:10,]
#
#

# ESCI functions ----------------------------------------------------------

# Conversion formula's from [Friedman(1982) and Wolf(1986)](http://www.soph.uab.edu/Statgenetics/People/MBeasley/Courses/EffectSizeConversion.pdf) and [Borenstein, Hedges & Rothstein (2009)](http://www.meta-analysis.com/downloads/Meta-analysis%20Converting%20among%20effect%20sizes.pdf)
f_d <- function(f,df1,df2){ifelse(df1>1,{2*sqrt(df1*f/df2)},{2*sqrt(f/df2)})}
f_r <- function(f,df1,df2){ifelse(df1>1,{sqrt((df1*f)/((df1*f) + df2))},{sqrt(f/(f + df2))})}
t_d <- function(t,df,n1=1,n2=1){((n1+n2)*t)/(sqrt(df)*sqrt(n1*n2))}
t_r <- function(t,df){sqrt(t^2/(t^2 + df))}
z_d <- function(z,N){(2*z)/sqrt(N)}
z_r <- function(z,N=NULL){# N>0 calculates effect size r, otherwise inverse Fisher Z-transform.
    ifelse(N>0,{return(cbind(rho=sqrt(z^2/(z^2 + N))))},{return(cbind(FisherZ.inv=tanh(z)))})}
X_d <- function(X,N,df=1){ifelse(df>1,{2*sqrt(X/N)},{2*sqrt(X/(N-X))})}
X_r <- function(X,N,df=1){ifelse(df>1,{sqrt(X/(X+N))},{sqrt(X/N)})}
r_z <- function(r,N=NULL){# N>=3 calculates effect size Z (Cohen's d, the standardised mean), otherwise Fisher Z-transform.
    ifelse(N>=3,{return(cbind(Standardised.Mean=(atanh(r)/(1/sqrt(N-3)))))},{return(cbind(FisherZ=atanh(r)))})}
r_d <- function(r){sqrt((4*(r^2))/(1-r^2))}
d_r <- function(d,n1=NULL,n2=NULL){
    #Borenstein, Hedges & Rothstein (2009):
    #The correction factor (a) depends on the ratio of n1 to n2, rather than the absolute values of these numbers.
    #Therefore, if n1 and n2 are not known precisely,use n1 = n2, which will yield 4.
    a <- 4
    if(!all(is.null(n1),is.null(n2))){a <- (n1+n2)^2/(n1*n2)}
    return(sqrt(d^2/((a+d^2))))
}
t_pr <- function(t,N,K=NULL){# partial correlation from t, see http://www3.nd.edu/~rwilliam/stats2/l02.pdf
    return(t/sqrt(t^2-(N-K-1)))}
b_t     <- function(b,SDb=NULL){b/SDb}
logOR_d <- function(logOR){logOR*(sqrt(3)/pi)}
d_logOR <- function(d){d*(pi/sqrt(3))}

cor.test.fisherZ <- function(r1=NULL,r2=NULL,n1=NULL,n2=NULL,p=TRUE){
    if((dim(as.matrix(r1))[2]==2)&(dim(as.matrix(r2))[2]==2)){
         r1 <- cor(r1[,1],r1[,2],use="pairwise.complete.obs")
         r2 <- cor(r2[,1],r2[,2],use="pairwise.complete.obs")
    } else{
     if(all((dim(as.matrix(r1))!=1))&all(dim(as.matrix(r2))!=1)){
         disp(message = "r1 and r2 each need to be:", header = "cor.test.fisherZ", footer = FALSE)
         disp(message = "- Either a single numerical value representing a correlation,", header = FALSE, footer = FALSE)
         disp(message = "- Or a 2 column matrix from which a correlation r1 and r2 can be calculated", header = FALSE)
     }
    }
    z <- ((atanh(r1)-atanh(r2))/((1/(n1-3))+(1/(n2-3)))^0.5)
    if(p){p<-2*(1-pnorm(abs(z)))} else {p=NULL}
    stat.test <- structure(list(statistic = z,
        method="Fisher r-to-Z transformed test for 2 independent correlations",
                                parameter=n1+n2,
                                p=p)
                           )
    class(stat.test) <- "htest"
    return(stat.test)
}


del2lam <- function(delta,N){# Change an observed t-value (delta) to a standardised mean (difference)
    if(length(N)>1){
        n1 <- N[1]
        n2 <- N[2]
    } else {
        n1 <- floor(N/2)
        n2 <- ceiling(N/2)
        if(N%%2>0){n1<-n1+1}
    }
    return(delta * sqrt((n1 * n2)/(n1 + n2)))
}

lam2del <- function(lambda,N){# Change an observed standardised mean (lambda) to a t-value (delta)
    if(length(N)>1){
        n1 <- N[1]
        n2 <- N[2]
    } else {
        n1 <- floor(N/2)
        n2 <- ceiling(N/2)
        if(N%%2>0){n1<-n1+1}
    }
    return(lambda * sqrt((n1 + n2)/(n1 * n2)))
}
#
# p2del <- function(p,N=NULL, df, prediction="two.sided"){
#   if(!is.null(N)){
#   if(length(N)>1){
#     df <- sum(N)-2
#   } else {
#     df <- N-1
#   }}
#   if(prediction=="two.sided") p<-p/2
#   return(delta = qt(p, df))
# }

sev.t <- function(t.obs,t.df,t.crit,discrepancy){
    require(MBESS)

    l  <- sapply(list(t.obs,discrepancy),length)
    if(l[2]!=1&l[2]<=10){
        cat("\nExpected a discrepancy vector of length == 1, or length >= 2\nWill use a 100 step sequence representing 95%CI for: t(",t.df,") = ",t.obs,"\n")
        CI <- conf.limits.nct(ncp=t.obs, df=t.df, conf.level=.95)
        discrepancy <- seq(1.5*CI$Lower.Limit,1.5*CI$Upper.Limit,length=100)
        l[2] <- 100
    }
    if(l[1]<l[2]){
        t.obs  <- rep(t.obs,length(discrepancy))
    }
    return(pt(t.obs,Ori.df,ncp=discrepancy,lower.tail=(t.crit<0)))
}


php.t = function(P, df, prediction="two.sided", alpha=.05) {# 'p-value based' post-hoc power
    # Lenth, R. (2007). Post hoc power: tables and commentary. http://www.stat.uiowa.edu/files/stat/techrep/tr378.pdf
    ifelse(prediction=="two.sided",two.tailed=TRUE,two.tailed=FALSE)

    if (two.tailed) {
        delta = qt(1 - P/2, df)
        cv = qt(1 - alpha/2, df)
        power = 1 - pt(cv, df, delta) + pt(-cv, df, delta)
    }
    else {
        delta = qt(1 - P, df)
        cv = qt(1 - alpha, df)
        power = 1 - pt(cv, df, delta)
    }
    power
}


posthocPOWer <- function(inference){

    if(inference$stat.type=="z"){
        posthocPOW <- if(inference$compare=="two.sided"){
            (1 - pnorm(abs(inference$stat.crit)) + pnorm(-1*abs(inference$stat.crit)))
        } else {
            (1 - pnorm(abs(inference$stat.crit)))
        }
    }

    if(inference$stat.type=="t"){
        posthocPOW <- if(inference$compare=="two.sided"){
            (1 - pt(abs(inference$stat.crit), inference$stat.df1,ncp=inference$stat.ncp) + pt(-1*abs(inference$stat.crit), inference$stat.df1,ncp=inference$stat.ncp))
        } else {
            (1 - pt(abs(inference$stat.crit), inference$stat.df1))
        }
    }
    if(inference$stat.type=="f"){
        posthocPOW  <- (1-pf(abs(inference$stat.crit), inference$stat.df1, inference$stat.df2))
    }
    if(stat.type=="chisq"){
        posthocPOW  <- (1 - pchisq(abs(inference$stat.crit), inference$stat.df1))
    }
    return(data.frame(posthocPOW=posthocPOW))
}

convertES <- function(inference, keepSign=TRUE){

    #   if(length(inference$stat.type)!=1){
    #     stop('\nChoose 1 distribution out of: ',paste(inference$stat.type,collapse="  "))
    #   } else {
    #     if(!any(tolower(c("z","r","pr","t","f","chisq","chi.sq","X^2","X2"))%in%inference$stat.type)){
    #       warning(paste("Unknown test statistic. Choose from: Z, t, F, chisq"),immediate.=T)
    #       ESconvert<-as.data.frame(t(rep(NA,6)))
    #       colnames(ESconvert) <- list("ES.d","ES.d.ciL","ES.d.ciU","ES.r","ES.r.ciL","ES.r.ciU")
    #       return(ESconvert)
    #     }
    #   }
    #
    #   if(any(is.null(inference$stat.ncp)||(inference$stat.ncp==0)||is.na(inference$stat.ncp))){
    #     ESconvert<-as.data.frame(t(rep(NA,6)))
    #     colnames(ESconvert) <- list("ES.d","ES.d.ciL","ES.d.ciU","ES.r","ES.r.ciL","ES.r.ciU")
    #     return(ESconvert)
    #   }
    #

    NCPvars <- list(ncp=inference$stat.ncp,ciL=inference$stat.ncp.ciL,ciU=inference$stat.ncp.ciU)

    ES.r    <- ldply(NCPvars,function(ncp) data.frame(rho=eval(parse(text=paste0('any2r(x=',ncp,',df1=',inference$stat.df1,',df2=',inference$stat.df2,',N=', (inference$stat.n1+inference$stat.n2),',esType="',inference$stat.type,'")'))) ))


    #    if(inference$stat.type=='z'){ CItext <- list(test='z',df=', (inference$stat.n1+inference$stat.n2)')}
    #   if(inference$stat.type=='t'){ CItext <- list(test='t',df=', inference$stat.df1')}
    #   if(inference$stat.type=='f'){ CItext <- list(test='f',df1=', inference$stat.df1',df2='inference$stat.df2')}
    #   if(inference$stat.type=='chisq'){ CItext <- list(test='X',df=',  (inference$stat.n1+inference$stat.n2), inference$stat.df1)')}
    #   if(inference$stat.type=='b'){ CItext <- list(test='b',df=', sum(inference$stat.N))')}
    #   if(inference$stat.type=='OR'){ CItext <- list(test='b',df=', sum(inference$stat.N))')}

    ES.d  <- ldply(ES.r$rho,function(ncp) data.frame(delta=r_d(ncp)))
    ES.OR <- ldply(ES.d$delta,function(ncp) data.frame(OR=d_logOR(ncp)))
    ES.fZ <- ldply(ES.r$rho,function(ncp) data.frame(fZ=atanh(ncp)))
    ES.v  <- ldply(ES.r$rho,function(ncp) data.frame(v=v.stat(n=(inference$stat.n1+inference$stat.n2),p=inference$stat.ncp.pval,Rsq=ncp^2)))

    ES.d$.id<-ES.r$.id
    ES.OR$.id<-ES.r$.id
    ES.fZ$.id<-ES.r$.id
    ES.v$.id<-ES.r$.id

    #   logORvars  <- list(ncp='ES.d$delta[[1]]',ciL='ES.d$delta[[2]]',ciU='ES.d$delta[[3]]')
    #   FisherZvars<- list(ncp='ES.r$rho[[1]]',ciL='ES.r$rho[[2]]',ciU='ES.r$rho[[3]]')
    #vVars     <-  list(ncp='ES.r$rho[[1]]',ciL='ES.r$rho[[2]]',ciU='ES.r$rho[[3]]')



    #   tmp<-list(ES.d,ES.d.ciL,ES.d.ciU,ES.r,ES.r.ciL,ES.r.ciU)
    #   tmp[sapply(tmp,length)==0]<-NA
    #   ESconvert<-as.data.frame(t(unlist(tmp)))
    #   if(keepSign){
    #     ESconvert <- ESconvert * as.numeric(c(1,1,1,sign(inference[,2:4])))
    #   } else {
    #     if(all(inference[,2:4]<=0)){
    #       ESconvert <- abs(ESconvert)
    #       message("\nkeepSign=FALSE\nCohen's d: All values were negative, removed sign...\n")
    #     } else {
    #       if((inference[,2:3]<=0)&(inference[,4]>0)){
    #         ESconvert     <- ESconvert - as.numeric(-1,-1,1,0,0,0)
    #         ESconvert[,3] <- 0
    #         message("\nkeepSign=FALSE\nCohen's d: Lower CI and NCP were negative, changed positive Upper CI to 0 and converted to 'right sided' values...\n")
    #       }
    #       if((inference[,2]<=0)&(inference[,3:4]>0)){
    #         ESconvert[,2] <- 0
    #         message("\nkeepSign=FALSE\nCohen's d: Upper CI and NCP were positive, changed negative Lower CI to 0...\n")
    #       }
    #     }
    #   }
    #   colnames(ESconvert) <- list("ES.d","ES.d.ciL","ES.d.ciU","ES.r","ES.r.ciL","ES.r.ciU")
    #
    options(warn = -1)
    return(data.frame(cbind(dcast(melt(ES.d), 1 ~ variable+.id)[2:4],dcast(melt(ES.r), 1 ~ variable+.id)[2:4],dcast(melt(ES.OR), 1 ~ variable+.id)[2:4],dcast(melt(ES.fZ), 1 ~ variable+.id)[2:4],dcast(melt(ES.v), 1 ~ variable+.id)[2:4])))
}

check.args <- function(stat.type, stat.ncp, stat.df1, stat.df2, stat.n1, stat.n2, alpha, CL, compare, tails){# Check input arguments

    stat.type<-tolower(stat.type)
    if(any(stat.type%in%tolower(c("r","pr")))){stat.type[stat.type%in%tolower(c("r","pr"))]<-"rho"}
    if(any(stat.type%in%tolower(c("B","b")))){stat.type[stat.type%in%tolower(c("B","b"))]<-"b"}
    if(stat.type%in%tolower(c("X^2","X2","chi.sq"))){stat.type<-"chisq"}
    #  if(stat.type%in%tolower(c("binom"))){stat.type<-"binom"}
    if(stat.type%in%tolower(c("z","Z","normal","norm"))){
        stat.type<-"z"
        cat("\nAssuming Mean = 0 and SD = 1/sqrt(N)")}

    if(is.na(tails)|(tails%in%"")){tails <- "two.sided"} else {
        if(tails==1){tails <- "one.sided"}
        if(tails==2){tails <- "two.sided"}
    }

    if(is.na(compare)|(compare%in%"")){compare <- "is"} else {
        if(compare=="<"){compare <- "lesser"}
        if(compare==">"){compare <- "greater"}
    }

    if(length(stat.n1)==0){stat.n1<-NA}
    if(length(stat.n2)==0){stat.n2<-NA}

    if(length(stat.type)!=1){
        stop('\nChoose 1 distribution out of: ',paste(stat.type,collapse="  "))
    } else {

        if(!any(tolower(c("z","t","f","chisq","r","pr","b","or"))%in%stat.type)){
            warning(paste("Unknown test statistic. Choose from: Z, t, F, chisq, r, B, binom"),immediate.=T)
            return(decide <- data.frame(stat.type="unknown",
                                        stat.ncp=stat.ncp,
                                        stat.df1=stat.df1,
                                        stat.df2=stat.df2,
                                        stat.n1=stat.n1,
                                        stat.n2=stat.n2,
                                        stat.ncp.ciL=NA,
                                        stat.ncp.ciU=NA,
                                        ci.type="unknown",
                                        stat.crit=NA,
                                        stat.ncp.p=NA,
                                        prediction=compare,
                                        tails = tails,
                                        stat.crit.p=alpha,
                                        H0=NA,
                                        H1=NA,
                                        decide="unknown"))
        }

        if( (all(is.na(stat.df1)&is.na(stat.df1)&stat.type=="f"))||(any(is.null(stat.ncp)||is.na(stat.ncp)||stat.ncp==0))){
            warning('Not enough information to evaluate the test statitsic.')
            return(decide <- data.frame(stat.type=stat.type,
                                        stat.ncp=NA,
                                        stat.df1=stat.df1,
                                        stat.df2=stat.df2,
                                        stat.n1=stat.n1,
                                        stat.n2=stat.n2,
                                        stat.ncp.ciL=NA,
                                        stat.ncp.ciU=NA,
                                        ci.type=NA,
                                        stat.crit=NA,
                                        stat.ncp.p=NA,
                                        prediction=compare,
                                        tails = tails,
                                        stat.crit.p=alpha,
                                        H0=TRUE,
                                        H1=FALSE,
                                        decide="Accept H0"))
        }
        #     if(length(prediction)>1){
        #       cat("\nArgument prediction > 1.\nAssuming undirected test...\n")
        #       prediction <- "is"
        #     } else {
        #       prediction <- tolower(prediction)
        #       if(!any(c("two.sided","greater","less")%in%prediction)){
        #         #stop("Unknown prediction, choose from: two.sided, greater, less")
        #         prediction<-"is"
        #       }
        #     }

        stat.type=enquote(stat.type)
        compare=enquote(compare)
        tails=enquote(tails)
        return(list(stat.type=stat.type, stat.ncp=stat.ncp, stat.df1=stat.df1, stat.df2=stat.df2, stat.n1=stat.n1, stat.n2=stat.n2, alpha=alpha, CL=CL, compare=compare, tails=tails))
    }
}

#    stat.type=ML2.primary[[s]]$stat.type
#   stat.ncp=ML2.primary[[s]]$stat.ncp
#     stat.df=ML2.primary[[s]]$stat.df
#     stat.N=ML2.primary[[s]]$stat.N
#     alpha=ML2.primary[[s]]$stat.info$stat.params$alpha
#     CL=ML2.primary[[s]]$stat.info$stat.params$conf.level
#     tails = 2#stat.[[s]]info$stat.params$tails,
#     compare=ML2.primary[[s]]$stat.info$stat.params$alternative

get.inference <- function(vlist){
    #vlist<-check.args()

    eval(parse(text=( paste0(names(vlist),"=",sep=paste(vlist)))))
    if(all(is.na(stat.df1),is.na(stat.df2),is.na(stat.n1),is.na(stat.n1))){return(vlist)}
    #stat.type=c("z","t","f","chisq","r"), stat.ncp, stat.df1, stat.df2, stat.N, alpha=0.05, CL=0.95, compare="is",tails=2
    if(stat.type!="f"){stat.df2 <-NA}

    ID    <- NULL
    alpha <- c(alpha,(1-(alpha)))
    H0=TRUE
    H1=FALSE

    if(stat.type=="z"){
        if(tails=="two.sided"){
            alpha=c(alpha[1]/2,(1-(alpha[1]/2)))
        } else {
            cat("\nDirectional Z-test... adjusted symmetric Confidence Level accordingly.")
        }
        CL <- diff(alpha)
        x.crit <- qnorm(alpha)
        stat.p <- pnorm(stat.ncp)
        MoE    <- qnorm(CL,lower.tail=F)*(1/sqrt(stat.n1+stat.n2))
        stat.ncp.ciL <- stat.ncp - MoE
        stat.ncp.ciU <- stat.ncp + MoE
        ci.type <- "symmetric"
    }

    if(stat.type=="t"){
        if(tails=="two.sided"){alpha=c(alpha[1]/2,(1-(alpha[1]/2)))
        } else {
            if(CL>(1-2*alpha)){
                CL <- 1-2*alpha
                cat("\nDirectional t-test... adjusted symmetric Confidence Level accordingly.")}
        }
        x.crit <- qt(alpha,stat.df1)
        stat.p <- pt(stat.ncp,stat.df1)
        CI <- conf.limits.nct(ncp=stat.ncp, df=stat.df1, conf.level=CL)
        stat.ncp.ciL <- CI$Lower.Limit
        stat.ncp.ciU <- CI$Upper.Limit
        ci.type <- "symmetric"
    }
    if(stat.type=="f"){
        if(tails=="two.sided"){tails="one.sided"}
        if(!any(is.na(stat.df1),is.na(stat.df2))){
            x.crit <- qf(alpha,df1=stat.df1,df2=stat.df2)
            stat.p <- pf(stat.ncp,df1=stat.df1,df2=stat.df2)
            if(stat.ncp>x.crit[2]){
                CI      <- conf.limits.ncf(F.value = stat.ncp, conf.level = CL, df.1 = stat.df1, df.2 = stat.df2)
                ci.type <- "symmetric"
            } else {
                CI     <- conf.limits.ncf(F.value = stat.ncp, conf.level = NULL, alpha.lower=alpha[1], alpha.upper=0, df.1 = stat.df1, df.2 = stat.df2)
                #CI$Lower.Limit <- stat.ncp
                CI$Lower.Limit <- stat.ncp
                ci.type <- "asymmetric"
            }
            if(length(CI$Lower.Limit)==0){CI$Lower.Limit <- NA}
            stat.ncp.ciL <- CI$Lower.Limit
            stat.ncp.ciU <- CI$Upper.Limit
            compare <- "greater"
        } else {
            warning("The F distribution requires 2 degrees of freedom.")
        }
    }
    if(stat.type=="chisq"){
        if(tails=="two.sided"){tails="one.sided"}
        x.crit <- qchisq(alpha,df=stat.df1,ncp=0)
        stat.p <- pchisq(stat.ncp,df=stat.df1,ncp=0)
        if(abs(stat.ncp-(x.crit[2]+1))>1){
            CI     <- try.CATCH(conf.limits.nc.chisq(Chi.Square=stat.ncp, conf.level=CL, df=stat.df1))
            ifelse(length(CI)==2,{
                CI <- conf.limits.nc.chisq(Chi.Square=stat.ncp, conf.level=NULL, alpha.lower=alpha[1], alpha.upper=0, df=stat.df1)
                ci.type <- "asymmetric"},{
                    CI <- CI$value
                    ci.type <- "symmetric"}
            )
        } else {
            CI     <- conf.limits.nc.chisq(Chi.Square=stat.ncp, conf.level=NULL, alpha.lower=alpha[1], alpha.upper=0, df=stat.df1)
            #CI$Lower.Limit <- stat.ncp
            ci.type <- "asymmetric"
        }
        if(length(CI$Lower.Limit)==0){CI$Lower.Limit <- NA}
        stat.ncp.ciL  <- CI$Lower.Limit
        stat.ncp.ciU  <- CI$Upper.Limit
        compare <- "greater"
    }

    if(compare=="is"){
        H0       <- (stat.ncp >=  x.crit[1] & stat.ncp <=  x.crit[2])
        H1       <- (stat.ncp <   x.crit[1] | stat.ncp >   x.crit[2])
        ifelse(stat.ncp < 0,{ID<-1},{ID<-2})
    }
    if(compare=="greater"){
        H0     <- (stat.ncp <= x.crit[2])
        H1     <- (stat.ncp >  x.crit[2])
        ID <- 2
    }
    if(compare=="lesser"){
        H0     <- (stat.ncp >= x.crit[1])
        H1     <- (stat.ncp <  x.crit[1])
        ID <- 1
    }

    ifelse(H0, say<-"Accept H0",say<-"Reject H0")
    ifelse(ID==2,{stat.pval=1-stat.p},{stat.pval=stat.p})

    decide <- tbl_df(data.frame(stat.type=stat.type,
                                stat.ncp=stat.ncp,
                                stat.df1=stat.df1,
                                stat.df2=stat.df2,
                                stat.n1=stat.n1,
                                stat.n2=stat.n2,
                                stat.ncp.ciL=stat.ncp.ciL,
                                stat.ncp.ciU=stat.ncp.ciU,
                                ci.type=ci.type,
                                stat.crit=x.crit,
                                compare=compare,
                                tails = tails,
                                stat.crit.p=alpha,
                                stat.ncp.p=stat.p,
                                stat.ncp.pval=stat.pval,
                                H0=H0,
                                H1=H1,
                                decide=say)
    )
    return(decide[ID, ])
}

v.stat <- function(n,p,Rsq){
    require(hypergeo)
    #From Daniel Lakens' blog: http://daniellakens.blogspot.nl/2014/11/evaluating-estimation-accuracy-with.html
    #Below, I'm vectorizing the function so that I can plot curves.
    #The rest is unchanged from the vstat function by Stober-Davis & Dana.
    #If you want to use R unbiased, remove the # before the Rsq adjustment calculation below

    #Rsq = Re(1-((n-2)/(n-p))*(1-Rsq)*hypergeo(1,1,(n-p+2)*.5,1-Rsq))
    if (Rsq<=0) {Rsq = .0001}
    r = ((p-1)*(1-Rsq))/((n-p)*Rsq)
    if(r<0|is.na(r)){r = abs(r)}
    g = min(r,1,na.rm=T)
    if (g<.5001 && g>.4999) {g = .5001}
    z = (g - sqrt(g-g^2))/(2*g - 1)
    alpha = acos((1-z)/sqrt(1-2*z*(1-z)))
    v = Re((((2*cos(alpha)*gamma((p+2)/2))/(sqrt(pi)*gamma((p+1)/2)))*(hypergeo(.5,(1-p)/2, 3/2, cos(alpha)^2) - sin(alpha)^(p-1))))
    return(v)
}

sev.data <- function(RPPdata,study){
    SEV.ori<- with(RPPdata, sev.info(stat.type=stat.ori.type[[study]],
                                     stat.ncp=as.numeric(stat.ori.ncp[[study]]),
                                     stat.df=c(as.numeric(stat.ori.df1[[study]]), as.numeric(stat.ori.df2[[study]])),
                                     stat.N=as.numeric(stat.ori.N[[study]]), prediction=prediction.ori[[study]],
                                     mus=list(SEV.ori.rep=as.numeric(stat.rep.ncp[[study]])),
                                     compare=list(SEV.comp.df=c(as.numeric(stat.rep.df1[[study]]),as.numeric(stat.rep.df2[[study]]))))
    )
    SEV.rep <- with(RPPdata, sev.info(stat.type=stat.rep.type[[study]],
                                      stat.ncp=as.numeric(stat.rep.ncp[[study]]),
                                      stat.df=c(as.numeric(stat.rep.df1[[study]]), as.numeric(stat.rep.df2[[study]])),
                                      stat.N=as.numeric(stat.rep.N[[study]]),
                                      prediction=prediction.rep[[study]],
                                      mus=list(SEV.ori.rep=as.numeric(stat.ori.ncp[[study]])),
                                      compare=list(SEV.comp.df=c(as.numeric(stat.ori.df1[[study]]),as.numeric(stat.ori.df2[[study]]))))
    )

    return(list(SEV.ori,SEV.rep))
}

sev.info <- function(stat.type=c("Z","t","F","chisq"),stat.ncp, stat.df, stat.N ,CL=.95, alpha=.05, prediction=c("two.sided","greater","less"),mus=list(SEV.0=0),compare=list(comp.df=1000)){
    require(MBESS)
    #   stat.ncp <- as.numeric(stat.ncp)
    #   stat.df[[1]]  <- as.numeric(stat.df[1])
    #   stat.df[[2]] <- as.numeric(stat.df[2])
    #   stat.N <- as.numeric(stat.N)
    #   CL       <- as.numeric(CL)
    #   alpha    <- as.numeric(alpha)

    stat.type<-tolower(stat.type)
    if(stat.type%in%tolower(c("X^2","X2","chi.sq"))){stat.type<-"chisq"}

    if(length(stat.type)!=1){
        stop('\nChoose 1 distribution out of: ',paste(stat.type,collapse="  "))
    } else {
        prediction <- tolower(prediction)
        if(!any(tolower(c("z","t","f","chisq","chi.sq","X^2","X2"))%in%stat.type)){
            stop("Unknown test statistic, choose from: Z, t, F, chisq")
        }

        #       if(length(prediction)>1){
        #         cat("\nArgument prediction > 1.\nAssuming undirected test...\n")
        #         prediction <- "two.sided"
        #       } else {
        #         if(!any(c("two.sided","greater","less")%in%prediction)){
        #           #stop("Unknown prediction, choose from: two.sided, greater, less")
        #           prediction <- "two.sided"
        #         }
        #       }

        if(stat.type!="f"){
            stat.df <- c(stat.df[1],NA)
        }

        infer      <- get.Stats(stat.type=stat.type,stat.ncp=stat.ncp,stat.df=stat.df,alpha=alpha,prediction=prediction)
        infer.comp <- get.Stats(stat.type=stat.type,stat.ncp=mus[[1]],stat.df=compare[[1]],alpha=alpha,prediction=prediction)

        if(is.na(infer$stat.ncp.ciL)){infer$stat.ncp.ciL<-0}
        if(is.na(infer$stat.ncp.ciU)){infer$stat.ncp.ciU<-0}
        x     <- seq(infer$stat.ncp.ciL,infer$stat.ncp.ciU,by=.01)
        x     <- sort(c(x,infer$stat.ncp.ciL,infer$stat.ncp.ciU,0))

        #     if(stat.type=="z"){
        #       CI    <- c(Lower.Limit=-1*infer$stat.crit,Upper.Limit=infer$stat.crit)
        #       x     <- seq(CI$Lower.Limit*1.5,CI$Upper.Limit*1.5,by=.01)
        #       x     <- sort(c(x,stat.ncp,CI$Lower.Limit,CI$Upper.Limit))
        #       y.ncl <- dt(x, dnorm=(ncp=CI$Lower.Limit)
        #       y.ncu <- dt(x, dnorm=(stat.df[1]),ncp=CI$Upper.Limit)
        #       y.ncp <- dt(x, df=(stat.df[1]),ncp=stat.ncp)
        #       POW   <- ifelse(prediction=="two.sided",{(1 - pt(abs(infer$stat.crit), stat.df[1], x) + pt(-1*abs(infer$stat.crit), stat.df[1], x))},{(1 - pt(abs(infer$stat.crit), stat.df[1], x))})
        #       SEV.crit <- pt(stat.ncp,stat.df[1],ncp=infer$stat.crit,lower.tail=(infer$stat.crit<0))
        #       SEV.obs  <- pt(stat.ncp,stat.df[1],ncp=stat.ncp, lower.tail=(infer$stat.crit<0))
        #     }

        if(stat.type=="t"){
            x.d   <- t_d(x,stat.df[1])
            x.r   <- t_r(x,stat.df[1])
            y.ncl <- dt(x, df=(stat.df[1]),ncp=infer$stat.ncp.ciL)
            y.ncu <- dt(x, df=(stat.df[1]),ncp=infer$stat.ncp.ciU)
            y.ncp <- dt(x, df=(stat.df[1]),ncp=stat.ncp)
            y.SEV.obs  <- pt(rep(stat.ncp,length(x)),stat.df[1],ncp=x,lower.tail=infer$H1)
            y.POW.post <- if(prediction=="two.sided"){
                1 - pt(rep(abs(infer$stat.crit),length(x)), stat.df[1], ncp=x) + pt(rep(-1*abs(infer$stat.crit),length(x)), stat.df[1], ncp=x)
            } else {
                (1 - pt(rep(abs(infer$stat.crit),length(x)), stat.df[1], ncp=x))
            }
            POW.post <- if(prediction=="two.sided"){
                (1 - pt(abs(infer$stat.crit), stat.df[1],ncp=stat.ncp) + pt(-1*abs(infer$stat.crit), stat.df[1],ncp=stat.ncp))
            } else {
                (1 - pt(abs(infer$stat.crit), stat.df[1]))
            }

            SEV.x   <- c(SEV.obs=stat.ncp,SEV.ciL=infer$stat.ncp.ciL,SEV.ciU=infer$stat.ncp.ciU,SEV.crit=infer$stat.crit,SEV.50= x[which.min(y.SEV.obs>=.5)],unlist(mus))
            SEV.d   <- t_d(SEV.x, stat.df[1])
            SEV.r   <- t_r(SEV.x, stat.df[1])
            SEV.y   <- sapply(SEV.x, function(mu) pt(stat.ncp, stat.df[1], ncp=mu, lower.tail=infer$H1))
            if(length(compare)==1){
                SEV.comp.x  <- qt(p=pt(SEV.x,stat.df[1],lower.tail=infer$H1), compare[[1]][1], lower.tail=infer.comp$H1)
                SEV.comp.y  <- sapply(SEV.comp.x,function(mu) pt(stat.ncp, compare[[1]][1], ncp=mu, lower.tail=infer.comp$H1))
            }
        }

        if(stat.type=="f"){
            prediction <- "greater"
            ifelse(length(stat.df)==2,{
                x.d   <- f_d(x,stat.df[1],stat.df[2])
                x.r   <- f_r(x,stat.df[1],stat.df[2])
                y.ncl  <- df(x, df1=(stat.df[1]), df2=(stat.df[2]), ncp=infer$stat.ncp.ciL)
                y.ncu  <- df(x, df1=(stat.df[1]), df2=(stat.df[2]), ncp=infer$stat.ncp.ciU)
                y.ncp  <- df(x, df1=(stat.df[1]), df2=(stat.df[2]), ncp=stat.ncp)
                y.SEV.obs  <- pf(rep(stat.ncp,length(x)),stat.df[1],stat.df[2],ncp=x,lower.tail=infer$H1)
                y.POW.post <- (1 - pf(rep(abs(infer$stat.crit),length(x)), stat.df[1], stat.df[2],ncp=x))
                POW.post   <- (1 - pf(abs(infer$stat.crit), stat.df[1], stat.df[2]))

                SEV.x   <- c(SEV.obs=stat.ncp,SEV.ciL=infer$stat.ncp.ciL,SEV.ciU=infer$stat.ncp.ciU,SEV.crit=infer$stat.crit,SEV.50= x[which.min(y.SEV.obs>=.5)],unlist(mus))
                SEV.d   <- f_d(SEV.x,stat.df[1],stat.df[2])
                SEV.r   <- f_r(SEV.x,stat.df[1],stat.df[2])
                SEV.y   <- sapply(SEV.x,function(mu) pf(stat.ncp,stat.df[1],stat.df[2],ncp=mu, lower.tail=infer$H1))

                if(length(compare)==1){
                    SEV.comp.x  <- qf(p=pf(SEV.x,stat.df[1],stat.df[2],lower.tail=infer$H1),df1=compare[[1]][1],df2=compare[[1]][2], lower.tail=infer.comp$H1)
                    SEV.comp.y  <- sapply(SEV.comp.x,function(mu) pf(stat.ncp,df1=compare[[1]][1],df2=compare[[1]][2],ncp=mu, lower.tail=infer.comp$H1))
                }
            },{
                stop("The F distribution requires 2 degrees of freedom.")
            })
        }

        if(stat.type=="chisq"){
            prediction <- "greater"
            x.d   <- X_d(x,stat.N,stat.df[1])
            x.r   <- X_r(x,stat.N,stat.df[1])
            y.ncl <- dchisq(x, df=(stat.df[1]),ncp=infer$stat.ncp.ciL)
            y.ncu <- dchisq(x, df=(stat.df[1]),ncp=infer$stat.ncp.ciU)
            y.ncp <- dchisq(x, df=(stat.df[1]),ncp=stat.ncp)
            y.SEV.obs  <- pchisq(rep(stat.ncp,length(x)),stat.df[1],ncp=x,lower.tail=infer$H1)

            y.POW.post <- (1 - pchisq(rep(abs(infer$stat.crit),length(x)), stat.df[1],ncp=x))
            POW.post   <- (1 - pchisq(abs(infer$stat.crit), stat.df[1], x))

            SEV.x   <- c(SEV.obs=stat.ncp,SEV.ciL=infer$stat.ncp.ciL,SEV.ciU=infer$stat.ncp.ciU,SEV.crit=infer$stat.crit,SEV.50= x[which.min(y.SEV.obs>=.5)],unlist(mus))
            SEV.d   <-  X_d(x,stat.N,stat.df[1])
            SEV.r   <-  X_r(x,stat.N,stat.df[1])
            SEV.y   <- sapply(SEV.x,function(mu) pchisq(stat.ncp,stat.df[1],ncp=mu, lower.tail=infer$H1))
            if(length(compare)==1){
                SEV.comp.x  <- qchisq(p=pchisq(SEV.x,stat.df[1],lower.tail=infer$H1),compare[[1]][1],lower.tail=infer.comp$H1)
                SEV.comp.y  <- sapply(SEV.comp.x,function(mu) pchisq(stat.ncp,compare[[1]][1],ncp=mu,  lower.tail=infer.comp$H1))
            }
        }

    }

    return(list(inference= data.frame(stat=stat.type,df1=stat.df[1],df2=stat.df[2],alpha=alpha,stat.CI=CL,infer,POW.post=POW.post),
                severity = data.frame(cbind(SEV.x,SEV.d,SEV.r,SEV.y,SEV.comp.x,SEV.comp.y)),
                curves   = data.frame(stat.x=x,stat.d=x.d,stat.r=x.r,y.ncl=y.ncl,y.ncu=y.ncu,y.ncp=y.ncp,y.POW.post=y.POW.post,y.SEV.obs=y.SEV.obs)))
}

plotReplication <- function(SEV.ori, SEV.rep, d.axis=c("stat","d","r"), studyname="Study 1", pl.sevlabels=TRUE, pl.sevlabels.comp=FALSE, pl.power=TRUE, pl.crit=TRUE, pl.labels=TRUE, pl.connect=FALSE, pl.ci=FALSE){
    require(ggplot2)

    if(length(d.axis)!=1){
        cat("\nArgument discrepancy axis not 1.\nAssuming original test statistic test...\n")
        d.axis <- "stat"
    } else {
        if(!any(c("stat","d","r")%in%d.axis)){
            stop("Unknown discrepancy axis, choose from: stat, d, r")
        }
    }

    ori.dat <- SEV.ori[["curves"]]
    ori.t   <- SEV.ori[["inference"]]
    ori.sev <- SEV.ori[["severity"]]

    rep.dat <- SEV.rep[["curves"]]
    rep.t   <- SEV.rep[["inference"]]
    rep.sev <- SEV.rep[["severity"]]

    switch(d.axis,
           "stat" = dID <- 1,
           "d"    = dID <- 2,
           "r"    = dID <- 3
    )


    df <- data.frame(rbind(ori=ori.dat,rep=rep.dat),study=c(rep(paste("Original:",ori.t$decide),times=nrow(ori.dat)),rep(paste("Replication:",rep.t$decide),nrow(rep.dat))),power=c(rep(paste("Original: POW =",round(ori.t$POW.post,digits=2)),times=nrow(ori.dat)),rep(paste("Replication: POW =",round(rep.t$POW.post,digits=2)),nrow(rep.dat))),row.names=NULL,stringsAsFactors=F)

    test <- c(paste0("Original d(x0) = ",round(ori.sev[1,dID],digits=2),"\nCrit. ",ori.t$stat,"(",ori.t$df1,if(!is.na(ori.t$df2)){paste0(", ",ori.t$df2)},") = ",round(ori.t$stat.crit,digits=2)),paste0("Replication d(x0) = ",round(rep.sev[1,dID],digits=2),"\nCrit. ",rep.t$stat,"(",rep.t$df1,if(!is.na(rep.t$df2)){paste0(", ",rep.t$df2)},") = ",round(rep.t$stat.crit,digits=2)))

    df.sevr <- data.frame(rbind(ori=ori.sev,rep=rep.sev),study=c(rep(paste("Original:",ori.t$decide),times=nrow(ori.sev)),rep(paste("Replication:",rep.t$decide),nrow(rep.sev))),yend=c(ori.sev$SEV.y,rep.sev$SEV.y),yend.comp=c(ori.sev$SEV.comp.y,rep.sev$SEV.comp.y) )
    xmin   <- min(df[,1])

    if(pl.sevlabels){
        df.sev <- df.sevr[c(1,6,7,12,4,10), ]
        mlt    <- grep("Original",df.sev$study,fixed=T)
        xmin   <- min(df[,1])-((max(df[,1])-min(df[,1]))/5)
        df.sev$xend <- c(xmin,(xmin+min(df[,1]))/2,xmin,(xmin+min(df[,1]))/2,xmin,(xmin+min(df[,1]))/2)[rank(df.sev$yend)]
        col <- 1:6
        col[ mlt] <- rep("grey50",length(mlt))
        col[-mlt] <- rep("black",length(mlt))
        #     df.sev$xend        <- df.sev[,1]
        #     df.sev$xend[ mlt]  <- xmin
        #     df.sev$xend[-mlt]  <- (xmin+min(df[,1]))/2
        df.sev$dfrom       <- c(test,rev(test),test)
        pre <- character(nrow(df.sev))
        ifelse(ori.t$H1,{pre[ mlt]<-"SEV(mu>"},{pre[ mlt]<-"SEV(mu<="})
        ifelse(rep.t$H1,{pre[-mlt]<-"SEV(mu>"},{pre[-mlt]<-"SEV(mu<="})

        #labels <- paste0(pre,round(df.sev[ ,dID],digits=2),")==",round(df.sev$SEV.y,digits=2))
        labels <- paste0(pre,round(df.sev[ ,dID],digits=2),")")
    }

    if(pl.sevlabels.comp){
        df.sev.comp <- df.sevr[c(6,12), ]
        mlt    <- grep("Original",df.sev.comp$study,fixed=T)
        dst    <- abs(mean(c(ori.sev[,1]-rep.sev[,1])))*1.2
        df.sev.comp$xend        <- df.sev.comp$SEV.comp.x
        df.sev.comp$xend[mlt]   <- df.sev.comp$xend[mlt]+(-1.5*dst)
        df.sev.comp$xend[-mlt]  <- df.sev.comp$xend[-mlt]+(1.5*dst)
        df.sev.comp$dfrom       <- rev(test)
        df.sev.comp$disc <- df.sev.comp$SEV.comp.x
    }

    if(pl.crit){
        df.crit<- df.sevr[c(4,10), ]
        df.crit$dfrom <- test
        df.crit$disc<- df.crit[ ,1]
    }

    # Discrepance axis will be labelled according to d.axis
    df.sev$disc <- df.sev[ ,1]
    df$disc     <- df[ ,1]
    d.label <- paste0("Discrepancy in units of ",list(ori.t$stat,"Cohen's d","Effect Size r")[[dID]])
    verdict <- ifelse(rep.t$H0,paste0("ReplicationInference: mu <= mu[1]"),paste0("ReplicationInference: mu > mu[1]"))

    repP <-  ggplot(df,aes(group=study)) +
        geom_line(data=df,aes(x=disc,y=y.SEV.obs,color=study),size=3) +
        geom_point(data=df.sev,aes(x=disc, y=SEV.y,shape=dfrom,fill=dfrom),alpha=.5,size=7)

    if(pl.sevlabels){repP <- repP +
        annotate("text",  x = df.sev$xend+.05, y = df.sev$yend, label = labels, size=2, parse=T, colour=col)}
    #annotate("text",  x = c(xmin,(xmin+min(df[,1]))/2), y = c(1,1), label = c("Original","Replication"), cex=4,fontface=2)}
    #==",round(df.sev.comp$SEV.comp.y,digits=2)
    if(pl.sevlabels.comp){repP <- repP +
        annotate("text",  x = df.sev.comp$xend, y = df.sev.comp$SEV.comp.y, label = paste0("SEV(mu>",round(df.sev.comp$SEV.comp.x,digits=2),")"),parse=T,size=10) +
        geom_point(data=df.sev.comp,aes(x=disc, y=SEV.comp.y,shape=dfrom,fill=dfrom),alpha=.5,size=7) +
        geom_segment(data=df.sev.comp[1:4, ],aes(x=SEV.x, y=SEV.y, xend=SEV.comp.x, yend=SEV.comp.y),color="grey80")}

    if(pl.power){repP <- repP + geom_line(data=df,aes(x=disc,y=y.POW.post,group=power,colour=study,linetype=power))}

    if(pl.connect){repP <- repP + geom_segment(data=df.sev[1:4,],aes(x=disc, y=SEV.y, xend=rev(disc), yend=rev(SEV.y)),color="blue",alpha=.5)}

    if(pl.crit){repP <- repP + geom_point(data=df.crit,aes(x=disc, y=SEV.y,shape=dfrom),color="red",fill="red",alpha=.6,size=3)}

    if(pl.ci){repP <- repP + geom_point(data=df.crit,aes(x=disc, y=SEV.y,shape=dfrom),color="red",fill="red",alpha=.6,size=3)}

    if(pl.labels){repP <- repP +
        ggtitle(paste(studyname)) +
        annotate("text",x=(max(df$disc)-xmin)/3,xmin=xmin,xmax=max(df$disc), y = 1.1, label = verdict, size=5,parse=T) +
        ylab("POWER / SEVERITY") +
        xlab(paste0(d.label))}
    #seq(round(min(df$disc),digits=1),round(max(df$disc)
    repP <- repP +
        scale_x_continuous(breaks=round(c(min(df$disc),df.sev$disc,max(df$disc)),digits=2),labels= round(c(min(df$disc),df.sev[ ,dID],max(df$disc)),digits=2),limits=c(xmin,max(df$disc))) +
        scale_y_continuous(breaks=c(0,df.sev$SEV.y,1),labels= round(c(0,df.sev$SEV.y,1),digits=2),limits=c(0,1.1)) +
        scale_color_manual(values=c("grey","grey30"),guide=guide_legend("Study: N-P decision")) +
        scale_shape_manual(values=c(21:25),guide=guide_legend(expression(paste("Evaluate ",d(x[0]),":",mu[1] == (mu[0] + gamma))))) +
        scale_linetype_manual(values=c(2,2),guide=guide_legend("Post-hoc power")) +
        scale_size_manual(values=c(1,2),guide=F) +
        scale_fill_manual(values=c("grey","black"),guide=guide_legend(expression(paste("Evaluate ",d(x[0]),":",mu[1] == (mu[0] + gamma))))) +
        theme_bw(base_size = 12, base_family = "")

    return(repP)
}

any2r <- function(x, df1, df2, N, esType){# by C. Hartgerink
    esComp <- ifelse(esType=="t",
                     sqrt((x^2*(1 / df1)) / (((x^2*1) / df1) + 1)),
                     ifelse(
                         esType=="f",
                         sqrt((x*(df2 / df1)) / (((x*df2) / df1) + 1))*sqrt(1/df2),
                         ifelse(
                             esType=="r",
                             x,
                             ifelse(
                                 esType=="chisq",
                                 sqrt(x/N),
                                 ifelse(
                                     esType == "z",
                                     tanh(x * sqrt(1/(N-3))),
                                     NA
                                 )
                             )
                         )
                     ))
    return(esComp)
}


get.ESCI <- function(data, CL=.95){
    # Function to get ESCIs based on the test statistic
    require(MBESS)
    if(is.na(data$Value)){
        res <- rep(NA,times=4)
        dr  <- as.data.frame(rbind(rep(NA, times=6))) # Need N!
        names(dr) <- c("d","d.lo","d.hi","r","r.lo","r.hi")}
    else {
        if(data$Statistic == "F"){
            res <- conf.limits.ncf(F.value=data$Value, conf.level=CL, df.1=data$df1, df.2=data$df2)
            dr  <- cbind(d=f_d(data$Value, df1=data$df1, df2=data$df2), d.lo=f_d(res[[1]], df1=data$df1, df2=data$df2), d.hi=f_d(res[[3]], df1=data$df1, df2=data$df2), r=f_r(data$Value, df1=data$df1, df2=data$df2), r.lo=f_r(res[[1]], df1=data$df1, df2=data$df2), r.hi=f_r(res[[3]], df1=data$df1, df2=data$df2) )
        }
        if(data$Statistic == "t"){
            res <- conf.limits.nct(t.value=data$Value, conf.level=CL, df=data$df1)
            dr  <- cbind(d=t_d(data$Value, df=data$df1), d.lo=t_d(res[[1]], df=data$df1), d.hi=t_d(res[[3]], df=data$df1), r=t_r(data$Value, df=data$df1), r.lo=t_r(res[[1]], df=data$df1), r.hi=t_r(res[[3]], df=data$df1) )}

        if(data$Statistic == "Chi2"){
            res <- conf.limits.nc.chisq(Chi.Square=data$Value, conf.level=CL, df=data$df1)
            dr  <- as.data.frame(rbind(rep(NA, times=6))) # Need N!
            names(dr) <- c("d","d.lo","d.hi","r","r.lo","r.hi")
        }
        if(data$Statistic == "Z"){
            if(data$df1>=100){
                res <- conf.limits.nct(t.value=data$Value, conf.level=CL, df=data$df1)
            }
            dr  <- as.data.frame(rbind(rep(NA, times=6))) # Need N!
            names(dr) <- c("d","d.lo","d.hi","r","r.lo","r.hi")
        }
        if(data$Statistic == "r"){
            res <- ci.R(R=data$Value, conf.level=CL, N=(data$df1+2), K=1)
            dr  <- cbind(d=r_d(data$Value), d.lo=r_d(res[[1]]), d.hi=r_d(res[[3]]), r=data$Value, r.lo=res[[1]], r.hi=res[[3]] )}
    }
    statci <- as.data.frame(cbind(levels(data$Statistic)[data$Statistic],CL,data$Value,rbind(res[c(1,3)])))
    names(statci) <- c("stat","CL","stat.v","stat.lo","stat.hi")
    return(cbind(statci,as.data.frame(dr)))
}


# scale.R -------------------------------------------------------------

# # Three uses:
# #
# # 1. scale.RANGE(x)             Scale x to data range: min(x.out)==0;      max(x.out)==1
# # 2. scale.RANGE(x,mn,mx)       Scale x to arg. range: min(x.out)==mn==0;  max(x.out)==mx==1
# # 3. scale.RANGE(x,mn,mx,lo,hi) Scale x to arg. range: min(x.out)==mn==lo; max(x.out)==mx==hi
# #
# # Examples:
# #
# # Works on numeric objects
# somenumbers <- cbind(c(-5,100,sqrt(2)),c(exp(1),0,-pi))
#
# scale.RANGE(somenumbers)
# scale.RANGE(somenumbers,mn=-100)
#
# # Values < mn will return < lo (default=0)
# # Values > mx will return > hi (default=1)
# scale.RANGE(somenumbers,mn=-1,mx=99)
#
# scale.RANGE(somenumbers,lo=-1,hi=1)
# scale.RANGE(somenumbers,mn=-10,mx=101,lo=-1,hi=4)

scale.R <- function(x,mn=min(x,na.rm=T),mx=max(x,na.rm=T),lo=0,hi=1){
    x <- as.data.frame(x)
    u <- x
    for(i in 1:dim(x)[2]){
        mn=min(x[,i],na.rm=T)
        mx=max(x[,i],na.rm=T)
        if(mn>=mx){warning("Minimum (mn) >= maximum (mx).")}
        if(lo>=hi){warning("Lowest scale value (lo) >= highest scale value (hi).")}
        ifelse(mn==mx,{u[,i]<-rep(mx,length(x[,i]))},{
            u[,i]<-(((x[i]-mn)*(hi-lo))/(mx-mn))+lo
            id<-complete.cases(u[,i])
            u[!id,i]<-0
        })
    }
    return(u)
}

Rmd2htmlWP <- function(infile, outfile, sup = T) {
    require(markdown)
    require(knitr)
    mdOpt <- markdownHTMLOptions(default = T)
    mdOpt <- mdOpt[mdOpt != "mathjax"]
    mdExt <- markdownExtensions()
    mdExt <- mdExt[mdExt != "latex_math"]
    if (sup == T) {
        mdExt <- mdExt[mdExt != "superscript"]
    }
    knit2html(input = infile, output = outfile, options = c(mdOpt), extensions = c(mdExt))
}

# MULTIPLOT FUNCTION ------------------------------------------------------------------------------------------------------------------
#
# [copied from http://www.cookbook-r.com/Graphs/Multiple_graphs_on_one_page_(ggplot2)/ ]
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols:   Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multi.PLOT <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
    require(grid)

    # Make a list from the ... arguments and plotlist
    plots <- c(list(...), plotlist)

    numPlots = length(plots)

    # If layout is NULL, then use 'cols' to determine layout
    if (is.null(layout)) {
        # Make the panel
        # ncol: Number of columns of plots
        # nrow: Number of rows needed, calculated from # of cols
        layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                         ncol = cols, nrow = ceiling(numPlots/cols))
    }

    if (numPlots==1) {
        print(plots[[1]])

    } else {
        # Set up the page
        grid.newpage()
        pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

        # Make each plot, in the correct location
        for (i in 1:numPlots) {
            # Get the i,j matrix positions of the regions that contain this subplot
            matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

            print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                            layout.pos.col = matchidx$col))
        }
    }
}

# TRY … CATCH -------------------------------------------------------------------------------------------------------------------------

##================================================================##
###  In longer simulations, aka computer experiments,            ###
###  you may want to                                             ###
###  1) catch all errors and warnings (and continue)             ###
###  2) store the error or warning messages                      ###
###                                                              ###
###  Here's a solution  (see R-help mailing list, Dec 9, 2010):  ###
##================================================================##

# Catch *and* save both errors and warnings, and in the case of
# a warning, also keep the computed result.
#
# @title tryCatch both warnings (with value) and errors
# @param expr an \R expression to evaluate
# @return a list with 'value' and 'warning', where
#   'value' may be an error caught.
# @author Martin Maechler;
# Copyright (C) 2010-2012  The R Core Team
#
try.CATCH <- function(expr){
    W <- NULL
    w.handler <- function(w){ # warning handler
        W <<- w
        invokeRestart("muffleWarning")
    }
    list(value = withCallingHandlers(tryCatch(expr, error = function(e) e),
                                     warning = w.handler),
         warning = W)
}


